Basile Sauvage

Stefanie Hahmann

Georges-Pierre Bonneau

Gershon Elber

Detail Preserving Deformation of B-spline Surfaces with Volume Constraint

Keywords: B-spline surfaces, constrained deformation, volume preserving, levelof-detail editing, multiresolution analysis

Geometric constraints have proved to be helpful for shape modeling. Moreover, they are efficient aids in controlling deformations and enhancing animation realism. The present paper addresses the deformation of B-spline surfaces while constraining the volume enclosed by the surface. Both uniform and non-uniform frameworks are considered. The use of level-of-detail (LoD) editing allows the preservation of fine details during coarse deformations of the shape. The key contribution of this paper is the computation of the volume with respect to the appropriate basis for LoD editing: the volume is expressed through all levels of resolution as a trilinear form and recursive formulas are developed to make the computation efficient. The volume constrained is maintained through a minimization process for which we develop closed solutions. Real-time deformations are reached thanks to sparse data structures and efficient algorithms.

Introduction

Level-of-detail (LoD) editing of free form curves and surfaces is now established as a valuable modeling tool [START_REF] Finkelstein | Multiresolution curves[END_REF]SDD95,SDS96]. It is an attractive application of multiresolution methods, because it allows modification of the overall shape of a geometric model at any scale while automatically preserving all fine details.

In contrast to classical control-point-based editing methods where complex detail preserving deformations need to manipulate a large set of control points, multiresolution methods can achieve the same effect by manipulating only a few control points of some low resolution representation.

Nonetheless, there are application areas, including CAGD and computer animation, where constrained deformations are required. Obviously, the constraints offer additional and finer control over the deformations applied to curves and surfaces. In the past many researchers have explored constrained deformations of free form curves and surfaces. Linear constraints such as position, tangency, orthogonality, and symmetry [BB89, Fow92, Gle92, WW92, FB93] are generally related to direct shape manipulation. These constraints offer the advantage of efficient processing, allowing for interactive manipulation of the free form geometry. Non-linear constraints that are commonly considered are the area [START_REF] Elber | Multiresolution curve editing with linear constraints[END_REF][START_REF] Hahmann | Area preserving deformation of multiresolution curves[END_REF], the volume [START_REF] Rappoport | Volume-preserving free-form solids[END_REF], second order differential constraints such as convexity [START_REF] Kaklis | Convexity-preserving interpolatory parametric splines of nonuniform polynomial degree[END_REF][START_REF] Plavnik | Surface design using global constraints on total curvature[END_REF], curve constraints [CW92, GL96, PGL + 02], and first and second order fairing functionals such as the bending energy of a thin plate [CG91, GC95, FRSW87, HS92, [START_REF] Bonneau | Variational design of rational bézier curves and surfaces[END_REF][START_REF] Hahmann | Shape improvement of surfaces[END_REF]. Satisfying non-linear constraints, however, requires intense computational effort, so that their use for interactive shape manipulation is generally very limited.

In the context of multiresolution, detail preserving editing with only linear and area constraints has been studied for planar curves. [START_REF] Elber | Multiresolution curve editing with linear constraints[END_REF] showed that the enclosed area can be viewed as a linear constraint. Incorporated into a multiresolution free form editing environment, the method allows the manipulation of non-uniform B-spline curves at different scales. In [START_REF] Hahmann | Area preserving deformation of multiresolution curves[END_REF], a wavelet-based multiresolution formula of the area constraint has been developed for uniform Bspline curves. Both methods are designed for real-time deformations with up to a few thousand control points. Yet, a direct generalization to three-dimensional surfaces is not obvious because the growth of the complexity in the volume computations would make it unusable for any interactive manipulation. The latter, however, is essential for any surface editing method. The volume constraint is important for many applications -for example, when one designs an airplane's fuselage that is assumed to hold a fixed volume. Volume constraint also belongs to Lasseter's principles of animations [START_REF] Lassiter | Principles of traditional animation applied to 3d computer animation[END_REF], which states that volume preservation enhances the realism when deforming characters, in computer animations. Volume preserving shape deformation has been considered for free form solids [START_REF] Rappoport | Volume-preserving free-form solids[END_REF], FFD [START_REF] Hirota | Fast volumepreserving free form deformation using multi-level optimization[END_REF], space deformation [START_REF] Wolfram Von Funck | Vector field based shape deformations[END_REF], meshes [START_REF] Lipman | Volume and shape preservation via moving frame manipulation[END_REF], and multiresolution meshes [START_REF] Botsch | Multiresolution surface representation based on displacement volumes[END_REF][START_REF] Basile Sauvage | Volume preservation of multiresolution meshes[END_REF]. The volume preservation of smooth spline surfaces in a multiresolution basis has not been considered before.

This paper presents a method for interactive detail-preserving editing of Bspline surfaces with volume constraint. Both uniform and non-uniform B-spline basis functions are addressed. It makes the following three main contributions:

• It generalizes [START_REF] Elber | Multiresolution curve editing with linear constraints[END_REF][START_REF] Hahmann | Area preserving deformation of multiresolution curves[END_REF] to three-dimensional B-spline surfaces with volume preservation. The appearance of the deformation is intuitively controlled by two settings: the scale and the extent of the deformation. The scale defines a threshold: if the details are finer than this threshold, they are preserved during the deformation process. The extent defines a portion of the surface that is involved in the deformation: beyond this extent the surface does not change.

• It develops volume formulas that are consistent with the LoD editing of B-spline surfaces. A so-called two-scale basis is used for non-uniform Bspline. A wavelet based multiresolution (MR) basis is used for uniform B-splines. First, the volume is expressed in trilinear form in the coefficients describing the surface with respect to these bases. Then, recursive formulas are developed to efficiently compute the trilinear forms.

• The volume computations and constraint solvings are presented as closedform solutions, which allow interactive editing. It is supported by a precise analysis of the algorithmic issues.

The paper is organized as follows. Section 2 describes our approach. In Section 3 we give the basics on B-spline surfaces and we specify the interaction with uniform and non-uniform MR representations. Section 4 presents the volume computations and explains how to preserve the volume during the deformation process. Finally, we present some results in Section 5 and conclude in Section 6.

In order to make this paper concise and chronological, both uniform and nonuniform frameworks are treated at the same time. Skipping the sections 3.3 and 4.4 at first reading would lead the reader through the non-uniform setting only.

Overview

This paper aims to provide a robust model for real-time MR volume preserving deformation of B-spline surfaces. The problem can be stated as follows. Given two pieces of information:

• An initial surface S enclosing a volume Θ ref ,

• The displacement of a point on the surface or a control point, we seek to define a deformation ∆ satisfying two constraints:

• The volume enclosed by the deformed surface S + ∆ equals Θ ref ,

• The displacement of the point, while two settings control the appearance of the deformation:

• The scale below which the details of the shape are preserved,

• The extent defines the portion of the surface involved in the deformation. This problem is generally under-constrained since there are only one scalar constraint and one vector constraint versus many degrees of freedom (the position of the control points). Among the infinitely many possible solutions, the proposed process tend to minimize the impact of the volume constraint: the key step solves the minimum change of control points.

The solution proposed in this paper has been validated by an interactive editing process illustrated in Figure 1. It states the essential issues. Steps i to iν (Figure 1) gather the pre-processing and the adjustment of the settings. Steps ν and νi constitute the core of the deformation process where efficiency is decisive: the backwards arrow νi → ν means that these two steps are looped during interactive editing. The computational effort due to the volume preserving is divided between the shaded steps and meets two requirements: The data storage limitation (step iii, Figure 1) and the interactivity (step νi, Figure 1). Let us now explain the function of each step:

Step i. The preliminary step consists of loading the surface and building the appropriate data structures.

Step ii. The scale of the deformation is set: the geometric information at finer scales (i.e., the details) will be preserved during the deformation process. Accordingly, the following steps depend on this specified scale. The way of specifying it depends on the surface model. As explained in Section 3, this is done by choosing either some knot sub-sequences (non-uniform B-splines) or an MR level (uniform B-splines).

Step iii. At this stage some coefficients that are used for the volume computation (during step νi) are pre-computed and stored. As shown in Figure 1, this step has to be processed again only if one desires a new deformation scale. Hence, interactivity is not required at this stage and preferably, as much of the whole computation effort as possible should be invested here. This step is discussed in detail in Section 4.

Step iv. In order to control the portion of the surface involved in the editing process, the user sets the extent of the deformation. Our model needs no assumption about the way this is chosen (see discussion Section 4.2).

Step v. The surface is manipulated. Regardless of the surface interaction medium (drag-and-drop, skinning, collision detection), the outcome can be formulated as positional constraints. Denote this interaction outcome as a "displacement": either a control point or a point on the surface is displaced (see Section 3.4).

Step vi. The volume-preserving deformation is computed. It respects the target deformation (step ν) while modifying the surface at the specified scale (step ii) and in the specified extent (Step iν) as well as preserving the volume. We justify and explain, in Section 4, the explicit solution as a constrained minimization problem.

Multiresolution deformation for B-spline surfaces

In this section, the framework for the B-spline surfaces' manipulation is presented. In Section 3.1, the basics are reviewed and the notations for patched surfaces are set up. It defines the fine or detailed surfaces. Since the main issue is to define detail-preserving deformations, two different models are proposed:

• A multiresolution (MR) analysis for uniform patches (Section 3.3). This is a uniform B-spline wavelet scheme [SDS96]. As explained in Section 3.3, it provides a unified basis for the surface and the deformation. Then the shape can be intuitively manipulated through an approximating control polyhedron. Moreover, the basis changes using linear filters are efficient. In the ensuing discussion, a superscript e denotes the coarse level of resolution of all the variables.

• A two-scale analysis for non-uniform patches (Section 3.2). Two basis (fine and coarse) are defined. This is a simpler but less convenient framework based on knot insertions and removals. This choice is motivated by two factors: non-uniform B-spline wavelet schemes are more complex and less efficient than uniform ones. The advantage of the non-uniform setting lies in its ability to handle more general shapes than the uniform setting, i.e. possibly C 1 discontinuities. In the ensuing discussion, all the variables relative to the coarse basis are denoted with a covering tilde (' ').

Eventually, we define how the surface is manipulated (Section 3.4). Following the application we propose to manipulate either the control points or a point on the surface.

In the following bold indices denote double indices, bold letters denote points or vectors in the embedding space R 3 , and covering arrows denote high-dimensional vectors.

Tensor-product B-spline surfaces

A surface S is given by a set {S q } q of tensor-product B-spline patches. Each patch S q is defined as an element of a functional space V by

S q (u, v) = (mu-1,mv-1) i=(0,0) p i ϕ i (u, v) , (u, v) ∈ [ 0, 1] 2 , ( 1
)
where i = (i u , i v ) are double indices and where ϕ iu,iv (u, v) = B iu (u)B iv (v) are the tensor-product B-spline functions of degree d × d defined over the knot sequences

U = {u 0 , u 1 , . . . , u mu+d } and V = {v 0 , v 1 , . . . , v mv+d }.
These basis functions are piecewise polynomial with global C d-1 continuity and they span V . With respect to this basis, the coefficients p i = (x i , y i , z i ) ∈ R 3 are called control points. Our model remains valid if the degree d differs following the directions u and v but for brevity we assume it does not.

Since we need the surface S to be a priori closed, manifold and orientable, we assume the patches are joined along their boundaries. In order to easily manage C 0 continuity between patches, we add the following boundary-interpolating conditions that fix (clamp) the so-called boundary-knots:

u 0 = u 1 = . . . = u d = 0 = v 0 = . . . = v d and u mu = u mu+1 = . . . = u mu+d = 1 = v mv = . . . = v mv+d .
This implies that the boundary curves only depend on the boundary control points (e.g., S q (0, v) is completely defined by (p 0,iv ) 0≤iv<mv ). Moreover, we assume that the degrees and the knot sequences of both patches correspond along their common boundary. Hence, the C 0 continuity between two neighboring patches is ensured by only equating the boundary control points.

Two-scale editing for non-uniform B-splines

The main objective is to define an arbitrary coarse deformation ∆ for some patch S q in order to compute a deformed surface S q + ∆. The coarser ∆ is, the larger is the scale of the deformation. "Two-scale editing" means that a "fine-scale" basis is used for defining the surface S q while a "coarse-scale" basis is used for defining the deformation ∆. The coarse B-spline basis is defined from the fine one by removing knots in the sequences U and V.

Following [START_REF] Elber | Multiresolution curve editing with linear constraints[END_REF], we define ∆ as an element of some approximation space V ⊂ V . It is a spline space based on knot sub-sequences U = { u 0 , . . . , u e mu+d } ⊂ U and V = { v 0 , . . . , v e mv+d } ⊂ V, such that: • The boundary knots are preserved. Thus, the boundary-interpolating property is preserved.

• Some interior knots are removed in order to enlarge the scale of the deformation.

Hence, V is spanned by a set of coarse B-splines φj of degree d × d. The deformation ∆ ∈ V is written as

∆(u, v) = ( e mu-1, e mv-1) j=(0,0) δ j ϕ j (u, v) with δ j ∈ R 3 . (2)
By selecting the knots to be removed, the user chooses the deformation scale (Step ii, Figure 1): the more knots s/he removes, the larger is the support of φj , and the coarser is the deformation. Note that • The coarse basis is used to define the deformation -i.e., ( δ j ) j are the unknowns, but

• The fine basis is required to apply the deformation to the surface -(i.e., to compute S q + ∆) since S q / ∈ V .

Therefore, ∆ is expressed in the fine basis as

∆(u, v) = i δ i ϕ i (u, v) with (0, 0) ≤ i ≤ (m u -1, m v -1),
where δ i ∈ R 3 are computed from the δ j , using knot-insertions [START_REF] Farin | Curves and surfaces for CAGD: a practical guide[END_REF].

Multiresolution editing for uniform B-splines

Let S q be a tensor-product patch as defined by Equation (1), with uniform knot sequences. The uniform framework is nothing but a particular case of the nonuniform one. Hence, the coarse deformation ∆ could be defined the same way. Uniform B-spline basis functions, however, are more amendable to wavelet-based MR schemes [START_REF] Finkelstein | Multiresolution curves[END_REF]. Therefore, we choose to define ∆ in a B-spline wavelet basis because this approach is both more efficient and more convenient:

• Translation-invariant and level-invariant MR filters make the basis changes easy and efficient.

• S q and ∆ lie in the same space. Hence, the deformed surface S q + ∆ is instantly computed in the MR basis.

• MR analysis provides a meaningful coarse approximation of the shape (namely the control polyhedron, see Figure 6 left), whose manipulation is convenient.

• A single parameter (the MR level) controls the scale of the deformation.

This section is not intended to give a complete presentation of MR surfaces. One-dimensional B-spline wavelets are detailed in [FS94, SDS96] and a tensorproduct example is shown in [SDS96].

Building an MR surface analysis with n levels requires the knot sequences to be uniform and the number of intervals to be multiples of 2 n as well:

u i+1 -u i = 1 pu2 n , d ≤ i < p u 2 n + d = m u , and v i+1 -v i = 1 pv2 n , d ≤ i < p v 2 n + d = m v ,
where p u and p v are arbitrary integer values.

The MR analysis is based on a sequence of nested approximation spaces

V 0 ⊂ V 1 ⊂ • • • ⊂ V n .
Each space V j is spanned by some scaling functions (ϕ j i ) i with respect to which the coefficients p j i ∈ R 3 are called control points. For each multiresolution level j three detail spaces W j-1 0 , W j-1 1 and W j-1 2 make the complement of V j-1 in V j :

V j = V j-1 ⊕ W j-1 0 ⊕ W j-1 1 ⊕ W j-1 2 .
The spaces are spanned by the wavelets ψ j 0,i , ψ j 1,i and ψ j 2,i . The coefficients d j 0,i , d j 1,i and d j 2,i are called detail coefficients. Notations and properties are collected in Table 1.

Space Dimension

Basis functions Coefficients

V j (p u 2 j + d) × (p v 2 j + d) ϕ j i (u, v) i p j i i = p j W j 0 p u 2 j × (p v 2 j + d) ψ j 0,i (u, v) i d j 0,i i = d j 0 W j 1 (p u 2 j + d) × p v 2 j ψ j 1,i (u, v) i d j 1,i i = d j 1 W j 2 p u 2 j × p v 2 j ψ j 2,i (u, v) i d j 2,i i = d j 2 Table 1: Multiresolution spaces.
Figure 2: Tensor-product filter bank. From left to right: decomposition process (using filters A and B); From right to left: reconstruction process (using filters P and Q).

The fine space V n = V contains the detailed patch S q . Its basis functions ϕ n i = ϕ i are the tensor-product B-spline functions and the coefficients p n i = p i are the B-spline control points. Given the fine control points p n i defining S q (u, v) = p n i ϕ n i (u, v), the patch can be written in the MR basis at any level e, 0 ≤ e ≤ n, as follows:

S q (u, v) = i p e i ϕ e i (u, v) + n-1 j=e i d j 0,i ψ j 0,i (u, v) + i d j 1,i ψ j 1,i (u, v) + i d j 2,i ψ j 2,i (u, v) . (3)
The coarse control points p e i encode an approximation of the shape (low frequency). They are the vertices of a coarse polyhedron controlling the surface (see Figure 6 left). The details d j k,i encode higher frequencies: the larger j, the higher the frequency and the finer the details.

The main advantage with regard to our objectives is the ability to manipulate the global shape at an arbitrary level e through the control points p e while preserving the finer details: e sets the scale of the deformation (Step ii, Figure 1) -i.e., the lower e is, the coarser the deformation. The deformation ∆ lies in the approximation space V e and is defined by

∆(u, v) = i δ e i ϕ e i (u, v). (4)
Hence, computing S q + ∆ amounts to adding p e + δ e , both of them being vectors of coefficients in R 3 for the scaling functions in V e (see equations (3) and ( 4)).

The MR coefficients in (3) are computed thanks to the filterbank algorithm [START_REF] Mallat | A wavelet Tour of Signal Processing[END_REF]. It is a linear, recursive and invertible process computing the change of basis illustrated in Figure 2. It is based on two decomposition filters A j and B j , and two reconstruction filters P j and Q j for each level j, which are applied in a tensor-product manner on the matrices of coefficients. Given these filters in the shape of sparse matrices, the decomposition step at level j consists of computing the control points and the details at level j -1 from the control points at level j:

p j-1 = A j p j (A j ) T , d j-1 0 = B j p j (A j ) T , d j-1 1 = A j p j (B j ) T , d j-1 2 = B j p j (B j ) T ;
the reverse process is called reconstruction step:

p j = P j p j-1 (P j ) T + Q j d j-1 0 (P j ) T + P j d j-1 1 (Q j ) T + Q j d j-1 2 (Q j ) T .
Note that the deformation (4) could be exactly computed in the non-uniform framework by repeatedly removing every second knot from U (once for each level of resolution from n to e). Then V = V e and equations (2) and (4) are equivalent. This shows that the non-uniform model is more general. In compensation the uniform model provides the convenient decomposition (3) and the ensuing control polyhedron.

Surface manipulation

The present method is independent of the surface manipulating tool. The deformation may be driven either by direct user interaction in a modeler (e.g., the mouse's displacement in a drag-and-drop operation), by skinning, or via collision detection, etc. Therefore, we propose to formalize it by fixing a displacement t for either a control-point or a point on the surface.

A simple surface manipulation scheme is the direct displacement of control point ī and is equivalent to fixing one coefficient of the deformation:

• set δ ī = t in the non-uniform framework;

• set δ e ī = t in the uniform framework. This simple scheme is especially convenient in the uniform framework because the surface may be manipulated through the control polyhedron whose vertices are {p e i }. Hence, δ e ī is meaningful: it corresponds to the displacement of the vertex ī.

Alternatively, displacing a point on the surface with parameter (ū, v), i.e. moving S(ū, v) to a new location S(ū, v) + ∆(ū, v), is equivalent to fixing ∆(ū, v) to a given value t (e.g. obtained through a click-and-drag operation):

• Constrain j δ j ϕ j (ū, v) = t in the non-uniform framework;

• Constrain j δ e j ϕ e j (ū, v) = t in the uniform framework. The advantage of this alternative is found in the direct manipulation of the surface. It is more intuitive in the non-uniform framework since there is no control polyhedron. Moreover it allows the enforcement of exact positional constraints on the surface.

Volume preserving deformation

This section details the key contributions of our method to maintaining the volume while preserving the details during the surface deformation process. The first contribution is formulas for volume preservation that are consistent with the LoD editing frameworks (see Section 3). The second contribution is efficient algorithms for these formulas and for constraining the surface during the deformation process.

In Section 4.1, the computation of the enclosed volume as a trilinear form is presented. In section 4.2, we define volume-preservation as a constrained minimization problem. In Sections 4.3 and 4.4, efficient solutions are proposed (for one single patch) and the algorithmic issues are discussed. Both surface models (non-uniform and uniform) and both displacements (either a control point or a point on the surface) are considered. Eventually, in Section 4.5, we extend the results to several patches with C 0 continuity.

Volume enclosed by a fine surface

The signed volume enclosed by a parametric surface may be written as an integral over the domain [START_REF] Gonzalez-Ochoa | Computing moments of objects enclosed by piecewise polynomial surfaces[END_REF][START_REF] Elber | Multiresolution curve editing with linear constraints[END_REF]. The (signed) volume enclosed by a patched surface S equals the sum over the patches:

Θ (S) = patches q Θ (S q ) , (5)
where

Θ(S q ) = 1 v=0 1 u=0 z ∂x ∂u ∂y ∂v - ∂x ∂v ∂y ∂u du dv, (6)
for each patch

S q (u, v) =   x(u, v) y(u, v) z(u, v)  
. By substituting Formula (1) into (6) the volume becomes a trilinear form with respect to the control points' coordinates:

Θ(S q ) = i,j,k x i y j z k [ 0,1] 2 ϕ k ∂ϕ i ∂u ∂ϕ j ∂v - ∂ϕ i ∂v ∂ϕ j ∂u du dv (7)
where the sum runs over (0, 0) ≤ i, j, k ≤ (m u -1, m v -1). The integrals cannot be computed on-the-fly (i.e., during step νi, Figure 1), for efficiency reasons; nor can they be pre-computed and stored because of the expected memory size (m 3 u m 3 v integrals per patch). Hence, we propose to separate the parameters u and v:

ϕ k ∂ϕ i ∂u ∂ϕ j ∂v - ∂ϕ i ∂v ∂ϕ j ∂u du dv = θ u jukuiu θ v ivkvjv -θ u iukuju θ v jvkviv , (8)
where

θ u ijk = 1 0 B i (u) B j (u) B ′ k (u) du, (9)
θ v ijk = 1 0 B i (v) B j (v) B ′ k (v) dv . (10)
It is now possible to pre-compute and store m 3 u values θ u ijk and m 3 v values θ v ijk during the initialization step i (Figure 1). Moreover, the B-spline functions have local support (containing d+1 intervals). Hence, these values are stored in sparse data structures with size O(m u d 2 ) and O(m v d 2 ). The computation through ( 9) and (10) may be either analytical or via a numerical approximation.

Problem statement

We now propose a mathematical statement of the problem sketched in Section 2: how to define a relevant deformation that fits the constraints and the settings? Since we have only two constraints, displacement and volume, we propose to solve a constrained minimization problem. We minimize the L 2 norm of the deformation's coefficients ( δi or δ e i ). The displacement constraint was discussed in Section 3.4. The volume constraint is Θ(S q + ∆) = Θ(S q ). The first challenge is to write this constraint as a function of the free variables ( δi or δ e i ), i.e., finding formulas that fit, respectively, the two-scale basis or the MR basis. The second challenge is to build an algorithm that satisfies both the memory and interactivity requirements. In other words, the pre-computation (Step iii, Figure 1) and the on-line processing (Step νi, Figure 1) must be carefully balanced.

Above, we mentioned two settings. Scale control was discussed in Sections 3.2 and 3.3. In order to control the extent of the deformation, the user chooses a set F of "free" control points that will be modified in order to preserve the volume. In our editing tool, F is a neighborhood of the displaced point whose size is specified by the user.

While the setting of this "active domain" that can undergo deformation is orthogonal to the work presented here, it is a crucial step in any MR/LOD editing environment. The specification of this active domain should be interactive while hiding the internal representation (i.e. parametric domain) as much as possible from the end user. One such possibility will allow the user to sketch a closed region on the surface, region that will be employed internally, in the parametric domain of the surface, to identify all control points (i.e. degrees of freedom) that affects this region.

Volume is a non-linear constraint. For efficiency purposes and in order to avoid the use of non-linear optimizations, linearization techniques will be used but without loosening mathematical exactness. Following [START_REF] Elber | Multiresolution curve editing with linear constraints[END_REF] the volume constraint is linearized by separating the deformation according to the x, y and z-axes. This makes the processing more efficient. In the following sections we detail only the x-axis deformation. Symmetric formulas are obtained for the y and z-axes.

An alternative linearization of the volume constraint is a first-order Taylor expansion, as proposed in [START_REF] Hahmann | Area preserving deformation of multiresolution curves[END_REF]. One single minimization is needed (instead of three in the present method) but the volume is just approximated. Both methods differ in some specific geometric configurations. Most of our tests, however, show very similar results.

Non-uniform B-splines 4.3.1 Volume constraint in the two-scale basis

In order to define ∆ such that the volume is preserved we regard the volume Θ(S q + ∆) as a function of the variables δ i . Since the function Θ(.) is trilinear, separating the deformation by axes makes it linear in each axis and hence much easier to manage. Consider the case δ i = ( δx i , 0, 0) and t = (tx, 0, 0). Other axes are deduced by symmetry. By substituting Formulas (1) and (2) into Equation (6), we get

Θ(S q + ∆) = Θ(S q ) + i Θ Y Z (i) δx i , (0, 0) ≤ ∀ i ≤ ( m u -1, m v -1), where Θ Y Z (i) = j,k y j z k ϕ k ∂ ϕ i ∂u ∂ϕ j ∂v - ∂ ϕ i ∂v ∂ϕ j ∂u du dv. (11)
Then, the volume constraint Θ(S q + ∆) = Θ(S q ) for the x-axis becomes

i Θ Y Z (i) δx i = 0, ( 12
)
where the sum is in fact limited to the specified extent i.e. it runs over i ∈ F.

Thanks to the linearity of the integral and the sum operators, one can use the pre-computed θ u and θ v to compute the coefficients Θ Y Z (i): Firstly compute

Θ Y Z (i) = j,k y j z k θ u jukuiu θ v ivkvjv -θ u iukuju θ v jvkviv (13) for (0, 0) ≤ i ≤ (m u -1, m v -1)
, from which the Θ Y Z (i) are then computed, using a knot-removal algorithm based on Boehm's recurrence. The details are given in Appendix A.

Solving the minimization problem

All the terms of the minimization problem are now known. One has to solve it for both manipulation models discussed in Section 3.4: Displacing a control point or displacing a point on the surface. In both cases a closed form of the resulting deformation { δx i } is derived.

Displacing a control point

When displacing the point ī, one fixes δx ī = tx. Excluding ī from the set F of free control points, the volume constraint (12) becomes:

Θ Y Z (ī) tx + i∈F Θ Y Z (i) δx i = 0.
Hence, the minimization is min

{ e δx i } i∈F | δx i | 2 subject to Θ Y Z (ī) tx + i∈F Θ Y Z (i) δx i = 0.
Using the Lagrange multiplier λ [Cia88], the minimization problem is equivalent to solving the unconstrained min-max problem of max

λ min { e δx i } i∈F | δx i | 2 + λ Θ Y Z (ī) tx + i∈F Θ Y Z (i) δx i ⇔ -→ ▽ i∈F | δx i | 2 + λ Θ Y Z (ī) tx + i∈F Θ Y Z (i) δx i = 0 ⇔    2 δx i + λ Θ Y Z (i) = 0, ∀i ∈ F, Θ Y Z (ī) tx + Θ Y Z (i) δx i = 0. Having j∈F | Θ Y Z (j)| 2 = 0, ( 14
)
the system has a full-rank and the unique solution is

λ = 2 tx Θ Y Z (ī) j∈F | Θ Y Z (j)| 2
, and

δx i = -tx Θ Y Z (ī) j∈F | Θ Y Z (j)| 2 Θ Y Z (i), ∀i ∈ F. ( 15
)
Equation ( 14) approaches zero in the singular case where the set of support functions, in F, all provide vanishing volumetric values to the total volume and hence can clearly have no affect on the final shape's volume, not to say preserve it.

Displacing a point on the surface

When displacing a point with parameter (ū, v), the minimization is subject to the volume (12) and position constraints (see Section 3.4):

min

{ e δx i } i∈F | δx i | 2 subject to Θ Y Z (i) δx i = 0 ϕ i (ū, v) δx i = tx .
Using, once again, the Lagrange multipliers (λ and µ), the problem is reduced to

max λ,µ min { e δx i } i∈F | δx i | 2 + λ i∈F Θ Y Z (i) δx i + µ i∈F ϕ i (ū, v) δx i -tx ⇔ -→ ▽ i∈F | δx i | 2 + λ i∈F Θ Y Z (i) δx i + µ i∈F ϕ i (ū, v) δx i -tx = 0 ⇔          2 δx i + λ Θ Y Z (i) + µ ϕ i (ū, v) = 0, ∀i ∈ F, Θ Y Z (i) δx i = 0, ϕ i (ū, v) δx i = tx. Let α = j∈F | ϕ j (ū, v)| 2 j∈F | Θ Y Z (j)| 2 -   j∈F ϕ j (ū, v) Θ Y Z (j)   2 .
The system has a unique full-rank solution when α = 0, a solution that is equal to

λ = 2 tx α j∈F ϕ j (ū, v) Θ Y Z (j) , µ = -2 tx α j∈F | Θ Y Z (j)| 2 ,
and

δx i = -λ 2 Θ Y Z (i) - µ 2 ϕ i (ū, v), ∀i ∈ F. (16)
α can vanish in three different case:

1. when ϕ j (ū, v) vanish for all j, 2. when Θ Y Z (j) vanish for all j, or 3. when ϕ j (ū, v), ∀j and Θ Y Z (j), ∀j, seen as two vectors of size |F|, are colinear vectors.

In case 1, clearly the selection of the functions in F cannot support the displacement constraint. Similarly, in case 2, the selection of the functions in F cannot support the volume constraint. In case 3, these two sets of support functions are linearly dependent and hence cannot simultaneously support the two linearly independent constraints of the displacement and volume.

Algorithmic issues

Here we list the successive computations done during interactive Step νi (Figure 1) and discuss the complexity. We go into details for the x-axis only. The complete algorithm successively treats the x, y and z coordinates. The costs are expressed with respect to the degree d of the surface, the total number m u m v of control points and the number #F of free control points.

• Computing {Θ Y Z (i)} costs O(m u m v d 2 ) using Formula (13). This low cost is due to the sparsity of θ u and θ v .

• Computing { Θ Y Z (i)} from {Θ Y Z (i)} costs O(m u m v d
) using the knotsremoval algorithm (see Appendix A). This is due to the efficiency of the knots-removal step.

• Computing { δx i } i∈F costs O(#F) through Formula (15) or (16). Having closed forms (15) and ( 16) solving the minimization problem is decisive here for efficiency.

• Computing {δx i } from { δx i } using the knots-insertion algorithm costs O(m u m v d).

• Adding the deformation to the surface

x i ← x i + δx i costs O(m u m v ).
Moreover, #F ≤ m u m v and we can assume that degree d is bounded (in most cases d ≤ 5 is sufficient). Hence, the global complexity of Step νi (Figure 1) is O(m u m v ). In other words, it is linear with respect to the number of control points. This is essential for interactivity. In a typical drag-and-drop operation every small displacement of the mouse (step ν, Figure 1) implies the computation of a new volume preserving deformation (Step νi, Figure 1). Three key points allow the linear cost during the interaction: sparse data structures, closed forms solution to the minimization problems, and the ability to conduct the pre-computations.

Pre-computing { ϕ i (ū, v)} i∈F is also needed when displacing a point on the surface. It costs O( m u d 2 + m v d 2 + #F), during Step iii (Figure 1):

• Computing { B iu (ū)} iu and { B iv (v)} iv with the one-dimensional De Boor algorithm costs O( m u d 2 + m v d 2 ). • Computing ϕ i (ū, v) = B iu (ū) B iv (v) for i ∈ F costs O(#F).

Uniform B-splines

Uniform B-splines are a less general framework than non-uniform ones, for surface design (e.g., sharp features are not possible). However, it enables the use of efficient MR schemes based on wavelets. This section shows that the MR basis is both efficient and convenient for defining and controlling the volume preserving deformation. The values needed for the volume constraint can be pre-computed through a recursive process using the reconstruction filters P and Q.

Volume constraint in the MR basis

In the uniform framework, ∆ and S q are written in the same MR basis. Rewriting Equations ( 7) and (8) with superscripts denoting the level of resolution (see Section 3.3), the volume in the fine basis of V n is given by Θ(S q ) = i,j,k

x n i y n j z n k θ n,u jukuiu θ n,v ivkvjv -θ n,u iukuju θ n,v jvkviv . (17)
Assume the surface is decomposed at some fixed level e. A similar formula exists: Θ(S q ) = i,j,k

x e i y e j z e k θ e,u jukuiu θ e,v ivkvjv -θ e,u iukuju θ e,v jvkviv ,

where x e , y e and z e are the coordinates of all the MR coefficients, i.e., of p e , d e , d e+1 , . . . , d n-1 . This formula is obtained by substituting the MR representation (3) of S q into Equation (7).

Figure 3: Recursive processing of the tensors. From left to right: θ n,u , θ n-1,u and θ n-2,u (or θ n,v , θ n-1,v and θ n-2,v ).

The sets θ e,u = (θ e,u iu,ju,ku ) iu,ju,ku and θ e,v = (θ e,v iv,jv,kv ) iv,jv,kv are rank 3 tensors, and can represented by three-dimensional arrays (see Figure 3). The fine level tensors θ n,u and θ n,v (see Figure 3 left) have been pre-computed through Formulas (9) and (10), and stored during initialization (Step i, Figure 1). Then, θ e,u and θ e,v are computed through a recursive process (from left to right). The transposed filter (P j ) T and (Q j ) T are applied in a tensor product manner:

• Following the 3 axes on the yellow (light) block,

• following the 2 axes on the red (mid grey) blocks,

• following the 1 axis on the blue (dark) blocks,

• and green blocks (on the right, down and below) are invariant.

Proving that this process makes Equations ( 17) and (18) identical is simple but tedious. An intuitive explanation considers the transition from Equation to Equation (18):

• x e , y e and z e are obtained from x n , y n and z n through the filterbank algorithm (using filters A j and B j ),

• θ e,u and θ e,v are obtained from θ n,u and θ n,v through this process (using filters (P j ) T and (Q j ) T ),

• and the filters are related by

A j B j = P j Q j -1 .
Computing θ e,u and θ e,v through this process is quite expensive but it is made only once during Step iii (Figure 1, pre-computation and storage).

Consider an x-axis deformation δ e i = (δx e i , 0, 0). By combining Equations ( 18) and (4), the volume constraint Θ(S q + ∆) = Θ(S q ) becomes i

δx e i Θ e Y Z (i) = 0, (19)
where for (0, 0)

≤ i ≤ (p u 2 e + d -1, p v 2 e + d -1) Θ e Y Z (i) = j,k
y e j z e k θ e,u jukuiu θ e,v ivkvjv -θ e,u iukuju θ e,v jvkviv .

(20)

Solving the minimization problem

The solutions in the present uniform framework are close to the previous nonuniform framework. Hence, we do not detail the calculations.

Displacing a control point

Displacing the point ī fixes δx e ī = tx. Excluding ī from F, the volume constraint (19) becomes:

Θ e Y Z (ī) tx + i∈F Θ e Y Z (i) δx e i = 0.
Hence, the minimization is min

{ e δx i } i∈F |δx e i | 2 subject to Θ e Y Z (ī) tx + i∈F Θ e Y Z (i) δx e i = 0,
and the solution is

δx i = -tx Θ e Y Z (ī) |Θ e Y Z (j)| 2 Θ e Y Z (i), ∀i ∈ F, (21)
where the sum runs over j ∈ F and it is assumed not to equal zero.

Displacing a point on the surface

Including the position constraint (Section 3.4) and the volume constraint (19), one gets min

{ e δx i } i∈F |δx e i | 2 subject to Θ e Y Z (i) δx e i = 0, ϕ e i (ū, v) δx i = tx. Having α = j∈F |ϕ e j (ū, v)| 2 j∈F |Θ e Y Z (j)| 2 -   j∈F ϕ e j (ū, v) Θ e Y Z (j)   2 ,
never vanish, the unique solution is ∀i ∈ F,

δx e i = tx α     j∈F |Θ e Y Z (j)| 2   ϕ e i (ū, v) -   j∈F ϕ e j (ū, v) Θ e Y Z (j)   Θ e Y Z (i)   . (22)

Algorithmic issues

The pre-processing (Step iii, Figure 1) includes the computation of θ e,u and θ e,v and of { ϕ i (ū, v)} i∈F (when displacing a point on the surface).

Interactive Step νi (Figure 1) treats the x, y and z-axes successively. Treating the x-axis involves:

• Computing {Θ e Y Z (i)} that costs O(m u m v #F) through Formula (20). • Computing {δx e i } i∈F that costs O(#F) through Formula ( 21) or ( 22). • Adding the deformation to the surface x e i ← x e i + δx e i that costs O(#F). Hence, the global cost is O(m u m v #F). However, it is not an optimal bound because the computation of {Θ e Y Z (i)} depends on the level e and on the sparsity of θ e,u and θ e,v . If e = n it costs O(#F) only.

Multi-patched surfaces

Modeling a closed surface with a single tensor-product patch restricts the topology. Let {S q } q be a set of patches modeling a closed surface. We extended the previous results while preserving C 0 continuity between the patches. For example, in the uniform setting one has to just modify the volume constraint (19) in the following way. If the control point i belongs to more than one patch, i.e., it is a boundary or a corner point, the coefficient Θ e Y Z (i) is replaced by the sum over the patches to which the points belongs: patches q Θ e Y Z (i q , q), where i q is the index of the control point with respect to the patch q, and Θ e Y Z (i q , q) is computed through Formula (20). The deformation δx e i is then applied on both sides of the join. The complexity of the resulting algorithm remains the same. Figures 5 and6 for instance involve 3 patches.

While the presented scheme handles only linear constraints that preserve C 0 continuity across surfaces, additional linear constraints could be added to preserve higher order continuity across the surfaces. The problem of posing the necessary constraints to preserve G 1 or C 1 (or even higher order) continuity is beyond this work and was addressed by many, including the authors [Pet95, HB00, HB03]. We foresee no major obstacle to overcome in supporting higher order continuity across surfaces, once these additional constraints are imposed, or any other additional linear constraints.

Results

Figure 4 shows the deformation of a unit cube (6 bi-cubic patches). The surface on the left illustrates the selection of the extent of the deformation (red/dark region) as a neighborhood of the displaced point across the joins. The point's parameter in the upper patch is (ū, v) = (0.7, 0.8) and the displacement vector is d = (0.2, 0.2, 0.9) T . The two views of the deformed surface (middle and right) illustrate the linearization by separating the 3 axis: the added volume caused by the bump displacement is balanced by a cavity in the opposite direction.

Conclusion and future work

A method for the precise preservation of the volume enclosed by tensor product B-spline surfaces during an MR editing process has been presented. Both uniform and non-uniform settings have been addressed. Volume formulas for two-scale non-uniform B-spline deformation and for MR uniform B-spline deformation have been developed. The portion of the surface involved in the deformation is easily selected. Moreover, the fine details of the shape can be preserved through an appropriate choice of basis: a two-scale basis in the non-uniform setting and an MR basis in the uniform setting. A computation of the volume with respect to these bases is proposed. The volume constraint is then linearized and included in a minimization process. A closed form of the solution allows interactive editing.

It is supported by a precise analysis of the algorithmic issues. The results show that the volume constraint is of interest in order to enhance the deformation's realism. This makes our method suitable for the animation of soft objects.

There are several aspects that in our opinion are worth further investigation. We presented a volume preserving deformation method for the non-uniform case that rests on a two-scale basis. It would be interesting to investigate the use of a non-uniform B-spline wavelet basis [START_REF] Kazinnik | Orthogonal decomposition of nonuniform bspline spaces using wavelets[END_REF][START_REF] Lyche | Spline wavelets of minimal support[END_REF].

Another direction for future work could be the extension of our method to G 1 continuous surfaces. It is well known that G 1 continuity for surfaces of arbitrary topology is a quite difficult problem. In fact, at a patch corner with more than four patches, the so-called "vertex compatibility problem" has to be solved [START_REF] Peters | Biquartic c1-surface splines over irregular meshes[END_REF][START_REF] Hahmann | Polynomial surfaces interpolating arbitrary triangulations[END_REF]. In a classical tensor-product setting, where four patches meet at a corner, the G 1 continuity can be expressed by linear equations with respect to the control points. In all the other cases this constraint might be non-linear.

Figure 1 :

 1 Figure 1: Volume preserving editing process. Expensive steps are shaded.

Figure 4 :

 4 Figure 4: Deformation of the unit cube. Each of the six faces is a 15 × 15 uniform bi-cubic patch. Left: initial surface with the extent (red/dark region) and the displacement (arrow). Middle and right: two views after volume preserving deformation. Wireframe shows the initial position.

Figure 5 :

 5 Figure 5: Deformation of a non-uniform B-spline surface. Left: Initial surface. Middle: The volume is free. Right: The volume is preserved.

Figure 5

 5 Figure5shows the deformation of a non-uniform B-spline surface. The initial surface (left) is deformed by pulling a point on the surface (on the back of the neck). When the volume is unconstrained (middle), the result lacks realism if we expect the knight to be made of some soft non-stretch material. Volume

Figure 6 :

 6 Figure 6: Deformation of a uniform B-spline surface. Frome left to right : control polyhedrons; initial surfaces; the volume is free; the volume is preserved.

Figure 6

 6 Figure 6 shows the deformation of uniform B-spline surfaces. The initial surfaces (left) are deformed by pushing a control point on the back of the neck. Results without volume preserving (middle) and with volume preserving (right) are compared. One can see that even a small deformation affects the appearance. The bottom row presents the same deformation with some fine details added to the surface. Since these details are basically encoded by the MR details, they are preserved during the deformation.

Figure 7

 7 Figure7illustrates both the detail and volume preservation for non-uniform surfaces made up of two 14 × 15 bi-cubic patches. Three initial surfaces (left) differ in their details only. Manipulation is processed by constraining only two points on the surface. The same deformations (b,c) and (d,e) are applied to all of them, without (b,d) and with (c,e) volume preservation. The preservation of the details is convincing while the global shape is interactively modified to enforce the volume constraint.

Figure 7 :

 7 Figure 7: Column (a) shows the original shapes. In columns (b) and (d) two deformations are applied to the object in column (a) without volume preservation. Columns (c) and (e) shows the same deformations but this time with volume preservation.

Algorithm 1

 1 Knot removal without extra memory allocation. Input: n = m u , vector u contains U and Θ(0. . . n -1) contains Θ Y Z (., i v ) Output: m = m u , u(0 . . . m -1) contains U and Θ(0 . . . m -1) contains Θ Y Z (., i v ). m ← d + 1 // boundary knots are preserved for l = d + 1 to n -1 do if u(l) is preserved then Θ(m) ← Θ(l) u(m) ← u(l) m ← m + 1 else // u(l) is removed for r = 0 to d -1 do // then md -1 ≤ m d -1 + r ≤ m -2 Θ(m-d-1+r) ← u(l) -u(m-d-1+r) u(l+r) -u(m-d-1+r) Θ(m-d-1+r)+ u(l+r+1) -u(l) u(l+r+1) -u(m-d+r) Θ(m-d+r) end for // special case r = d Θ(m-1) ← u(l) -u(m-1) u(l+d) -u(m-1) Θ(m-1) + Θ(l) end if end for

Acknowledgments

The work was partially supported by the AIM@SHAPE Network of Excellence (FP6 IST NoE 506766) and the IMAG project MéGA.

A Knot removal algorithm

The volume preserving deformation of non-uniform B-spline surfaces as proposed in Section 4.3 requires efficient computation of Θ Y Z (i) from Θ Y Z (i). We propose a recursive algorithm based on Boehm's recurrence [START_REF] Boehm | Inserting new knots into a bspline curve[END_REF][START_REF] Farin | Curves and surfaces for CAGD: a practical guide[END_REF] that relates between the B-spline basis functions before and after a knot is inserted or removed.

For the purpose of clarity, first assume V = V and one single knot u l is removed from U, i.e. U = U ∪ {u l }. Thus, the coarse B-splines B k defined over the knot sequence U are related to the fine B-splines B k defined over U by

and B k = B k otherwise. Thanks to the linearity of the integral and of the sums (Equations ( 11) and ( 13)), this relation can be transfered between Θ Y Z (i) and Θ Y Z (i). The integral in Equation ( 11) may be written as

Hence, comparing Equations ( 11) and ( 13) implies

Using formula (23) for all i v one can compute each column Θ Y Z (., i v ) (relative to the coarse basis) from the column Θ Y Z (., i v ) (relative to the fine basis).

Then assume V = V and U ⊂ U. To each column Θ Y Z (., i v ) separately one has to apply Formula (23) for each removed knot, i.e., each knot in U \ U. The iterative Algorithm 1 is an example where no extra memory is needed.

Finally, whatever V and U: i) remove the knots u ∈ U \ U as explained above and ii) remove symmetrically the knots v ∈ V \ V for each row Θ Y Z (i u , .).

Algorithm 1 costs O(m u d). Hence, computing all the Θ Y Z (i) from the Θ Y Z (i) costs O(m u m v d).


Michel Brune

Julien Bernu

Christine Guerlin

Samuel Deléglise

Clément Sayrin

Sébastien Gleyzes

Stefan Kuhr

Igor Dotsenko

Jean-Michel M Raimond

Serge Haroche

‡ S Kuhr

§ I Dotsenko

Process tomography of field damping and measurement of Fock state lifetimes by quantum non-demolition photon counting in a cavity

Keywords: numbers: 03, 65, Ta, 42, 50, Pq

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The goal of quantum process tomography is to determine experimentally the matrix elements of the superoperator describing the evolution of a quantum system's density matrix [1]. This information is acquired by preparing a set of test states and monitoring their subsequent evolution. The method has been applied so far to spins in NMR experiments [2], to solid state qubits [3], to vibrational states of atoms [4] and to quantum gate operations [5,6]. We describe here process tomography applied to the photon number distribution of a relaxing field stored in a high-Q superconducting cavity, in which Fock states are used as test states.

Cavity field relaxation is described by a rate equation which, restricted to the photon number distribution P (n, t), is [START_REF] Haroche | Exploring the Quantum[END_REF] dP (n, t) dt = n ′ K n,n ′ P (n ′ , t).

(1)

Quantum electrodynamics predicts K n,n = -κ[(1 + n b )n + n b (n + 1)], K n,n+1 = κ(1 + n b )(n + 1), K n,n-1 = κn b n, all the other coefficients being 0. Here, κ is the damping rate of the field energy and n b the mean number of blackbody photons at temperature T . The time constant -1/K n,n is the lifetime of the n-Fock state. In this Letter, we report a complete experimental determination of the K n,n ′ coefficients.

Our experiment relies on a QND procedure [START_REF] Brune | [END_REF] to count the number n of photons stored in a cavity. It is based on the measurement of cavity-field-induced light-shifts on Rydberg atoms crossing the cavity one by one. It projects * Electronic address: brune@lkb.ens.fr † Present address: ETH Zurich, CH-8093 Zurich Switzerland ‡ Present address: IOTA, 91127 Palaiseau Cedex § Present address: Johannes Gutenberg Universität, Institut für Physik, Staudingerweg 7, D-55128 Mainz, Germany the field onto Fock states with high fidelity. By detecting long sequences of QND probe atoms along single field realizations, we follow the field evolution and observe the jumps between Fock states due to cavity damping [9,10]. By analyzing a large ensemble of field trajectories, we now partially reconstruct the super-operator describing the field relaxation process in the cavity and measure the lifetimes of individual Fock states which scale as 1/n at zero temperature [11]. This study provides insights into the physics of these highly non-classical states whose production by random [9,10] or deterministic [12] processes has recently been demonstrated.

Our setup [9,13] We consider fields with a negligible probability of having n > 7. In the Bloch pseudo-spin representation, the atomic state at the exit of C points along one out of 8 directions equally distributed in the equatorial plane of the Bloch sphere, corresponding to values of n varying from 0 (axis Ox) to 7. After leaving C, the atoms are submitted to a second π/2 pulse in R 2 with an adjustable phase φ with respect to that of R 1 . The atoms are detected by the field ionization counter D discriminating the states e and g. Measuring the atomic energy after R 2 amounts to detecting the atomic spin at the exit of C in the direction at an angle φ with Ox in the equatorial plane of the Bloch sphere. On average, we detect one atom every 0.24 ms (∼ 500 atoms detected during T c ).

Figure 1b shows a 3D histogram of the transverse atomic spin (components σ x and σ y ) after interaction with a coherent field in C. Each point in the Bloch sphere equatorial plane is obtained by measuring the average value of spin projections, on a sample of 110 consecutive atoms crossing C in a 26 ms time interval, much shorter than T c . About 700 atoms are sent across C, out of which we extract ∼ 600 atomic samples of 110 consecutive atoms. The procedure is repeated 2000 times, yielding about 10 6 spin measurements. The histogram clearly shows that the direction of the atomic spin is quantized.

After a measurement indicating a spin pointing towards a peak of this histogram, the field is projected onto the corresponding Fock state. This is checked by correlating two independent successive samples of 110 atoms along a single field realization. The first pins down n and the second remeasures it. Figure 1c shows the histogram of second measurements after selection of n = 3. It exhibits a main n = 3 peak with an n = 2 satellite due to field relaxation during the 26 ms time delay between the two measurements. This satellite can be suppressed by post-selection. Figure 1d displays the histogram of the intermediate results in sequences of 3 measurements for which the first and the last yield n = 3. The single peak reveals that, at the intermediate measurement time, the field contains exactly 3 photons. We use such single photon number peaks to calibrate the phase-shifts Φ(n).

Although the above method allows us to prepare Fock states and to observe qualitatively their evolution, it lacks the resolution required for a precise time analyzis.

The interval between two measurements (26 ms) is longer than the lifetime of n = 7 (18 ms). We can however analyze the data in a more efficient way, making better use of our atomic detection rate. On a single field realization, each atom detected along direction φ provides one bit of information j (j = 0 for e and j = 1 for g). After detecting N atoms, our knowledge of the field is described by an inferred photon number probability distribution p i N (n) linked to the initial distribution P 0 (n) by Bayes law:

p i N (n) = P 0 (n)Π N (n)/Z
where Z is a normalization and Π N (n) is the product of N functions p(j, φ|n), each describing the information provided by one atomic detection: p(j, φ|n) = 1/2[1 + (-1) j (A + B cos[Φ(n) + φ])] [10,15]. For successive atoms, we use different values of φ (-1.74, -0.87, 0 and 0.54 rad) chosen so that p(j, φ|n) is nearly maximal for n =6, 7, 0 or 1, respectively. The values of A and B, ideally 0 and 1, become -0.1 and 0.7, respectively, because of experimental imperfections. For N ∼ 100, p i N (n) converges to a Dirac peak corresponding to the photon number given by the atomic spin analysis.

Let us call P N (n) the ensemble average p i N (n) over many realizations in which the field is initially described by P 0 (n). As the detection process is QND, we have P 0 (n) = P N (n) for any N . In other terms, P 0 (n) is a fixed point of the transform P 0 (n) → P 0 (n)Π N (n)/Z . This property allows us to determine P 0 (n) by iteration of this transform starting with any initial non-vanishing distribution, for instance the flat one P fl (n) = 1/8. This method can be applied for determining P (n, t) at any time t by selecting in each sequence the sample of N detected probe atoms starting at this time.

We first reconstruct in this way the evolution of P (n, t) for a coherent field injected by S at t = 0 and relaxing in C. A measurement sequence, involving about 2750 atoms detected in 650 ms, is repeated 2000 times. We reconstruct P (n, t) with the above procedure using N = 25 atoms and 20 iterations. At each time t, we start the iteration with P fl (n). The temporal resolution is ∼ 6 ms, much shorter than T c /7.

Figure 2a presents P (n, t) versus t for n = 0 to 7 and Fig. 2b the time evolution of the average photon number n = n nP (n, t). According to theory, n is exponentially damped towards an offset corresponding to the blackbody background. The experimental decay (solid line) is indistinguishable from an exponential fit (red line), which yields T c = 132 ms in agreement with the independently determined value of the cavity damping time. The offset yields n b = 0.06, close to the theoretical value (0.05) of the blackbody field at the cavity temperature, 0.8 K. The insets in Fig. 2b present snapshots of P (n, t) at three different times with the corresponding Poisson fits. These histograms show, as theoretically predicted [START_REF] Haroche | Exploring the Quantum[END_REF], that the field remains coherent under the effect of damping (at the limit where blackbody effects are negligible). The dotted line in Fig. 2a presents a numerical solution of the theoretical rate equation using the above determined values for T c , n b and n at t = 0. It is in excellent agreement with the data (solid lines). We go now one step further. By monitoring the decay of selected Fock states, we determine the K n,n ′ coefficients without any a priori assumption about their values. The same experimental data is processed in two steps. First, we analyse separately the 2000 realizations of the experiment in order to select individual Fock states. For each sequence, we compute, after each atom detection, the new inferred photon number distribution p i (n, t) according to Bayes law [10]. We start the analyzis of each sequence with the Poisson distribution determined above. Between atoms, we evolve the estimated p i (n, t) according to the theoretical rate equation. This method gives, at each time, the best estimate of the actual photon number distribution in each realization. Except around quantum jumps, p i (n, t) is generaly peaked at a single photon number value n 0 . Whenever p i (n 0 , t) > 0.7 we assume that, within a good approximation, the n 0 Fock state is present in C at this time, which we take as origin (t = 0) for subsequent analysis of this Fock state decay.

In a second step, we gather all atomic data following the selection of a given n 0 , obtaining thus ensembles of Fock state-selected field realizations. We apply to each ensemble the iterative analysis described above, reconstructing for each value of n 0 the subsequent P n0 (n, t) distributions. As in the case of a coherent state, we use N = 25 atoms and 20 iterations starting with a flat initial distribution. Let us stress that this reconstruction procedure does not rely on any theoretical assumption about the form of the relaxation process. We made use of our theoretical knowledge of the K n,n ′ coefficients only in the first step of the data processing, in order to optimize the selection of the initial Fock states. Figure 3 shows in solid lines the reconstructed P n0 (n, t) distributions versus time for n 0 = 0 to 7 (a to h). In each frame, P n0 (n 0 , t) is, as expected, maximum at t = 0, its value giving the fidelity of the Fock state selection procedure. The other P n0 (n, 0) values are small. At long time (400 ms) the most probable photon number is always n = 0, reflecting the irreversible evolution of the field toward the thermal background close to vacuum. For n 0 = 0 (Fig. 3a) P 0 (0, t) decreases slightly below 1, while P 0 (1, t) reaches an equilibrium value close to 0.06. This describes the thermalization of the initially empty cavity. For n 0 = 1 (Fig. 3b) we observe the exponential decay of P 1 (1, t), together with the increase of P 1 (0, t), which describes the damping of a single photon into vacuum [9]. For n 0 > 1 (Fig. 3c to h), P n0 (n 0 , t) decreases exponentially at a rate increasing with n 0 (damping of the initial Fock state). The P n0 (n, t) functions with n = n 0 -1, n 0 -2, ..., 1 exhibit bell-shaped variations. They peak successively, reflecting the cascade of the photon number from n 0 down to vacuum. In order to extract the damping coefficients, we fit the first 20 ms of these curves to a solution of Eq. ( 1), leaving as free parameters the K n,n ′ and the initial P n0 (n, 0) values. The procedure is iterative. We get a first approximation of K n,n ′ with n and n ′ ≤ 1 using the data of Fig. 3a andb. We then determine the K n,n ′ with increasing indices by including progressively in the fits the data of Fig. 3c to h, optimizing at each step the previously determined parameters.

The obtained -K n,n values, which represent the decay rates of the n-Fock states, are shown versus n in Fig. 4a (in units of κ). As expected, they vary linearly with n. The solid straight line corresponds to the theory for n b = 0.06, while the dotted line shows the expected variation of K n,n at T = 0 K. This constitutes the first measurement of Fock states lifetime for n > 1, exhibiting clearly the expected 1/n variation [16]. Moreover, the departure of the experimental points from the dotted line shows that our procedure is precise enough to be sensitive to the small effect of the residual n b = 0.06 photon blackbody field on the lifetime of Fock states. The non-diagonal K n,n ′ coefficients are shown (Fig. 4b) in a 3D plot, in logarithmic scale. The big and small bars near the diagonal correspond to the K n,n+1 and K n,n-1 coefficients, respectively. The latter, which represent the thermal rates of photon upward jumps, are predicted to vanish for n b = 0. The logarithmic scale is convenient to display together the K n,n+1 and K n,n-1 coefficients which differ by about one order of magnitude for n b = 0.06. All other non-diagonal coefficients are 0 within noise. Figure 4c shows for comparison the corresponding theoretical coefficients for n b = 0.06.

The dotted lines in Fig. 3 are the result of a numerical integration of Eq. ( 1) using the values of K n,n ′ and P n0 (n 0 , 0) determined by our fit. The excellent agreement with the experiment over the full 400 ms time interval shows the accuracy of our method. This study demonstrates the power of QND photon number measurements to investigate the quantum behavior of a field stored in a cavity. It directly probes our theoretical understanding of field relaxation and clearly illustrates the high sensitivity of large photon number states to decoherence, their lifetime being (at T = 0 K) inversely proportional to their energy. The method is limited here to probing the rate equation for the photon number probability distribution. In order to determine the full super-operator of decoherence, we plan to monitor the evolution of the non-diagonal elements of the field density operator in the Fock state basis, using a time resolved quantum state reconstruction procedure demonstrated in Ref. [START_REF] Deléglise | Nature, in press and[END_REF].

FIG. 1 :

 1 FIG. 1: (a) Scheme of experimental set-up. (b) Histogram of the transverse atomic pseudo-spin after interaction with a coherent field in C (initial mean photon number 4.4). Photon numbers associated to each peak are given. (c) Spin histogram after selection of the n = 3 Fock state. The n = 2 peak is due to photon loss between selection and measurement. (d) Same histogram as in (c) conditioned to a post-selection measurement excluding events in which a photon is lost.

FIG. 2 :

 2 FIG. 2: Relaxation of a coherent state. (a) Evolution of the photon number probabilities P (n, t) (n = 0 . . . 7 according to the colour code defined in the inset). Black dotted lines are theoretical. (b) Average photon number versus t (solid black line) and exponential fit (thin red line). Insets show the photon number distributions (red histogram) and their Poisson fit (blue lines) at the three times shown by arrows.

FIG. 3 :

 3 FIG.3: Relaxation of Fock states. (a) to (h) Evolution of the photon number distributions Pn 0 (n, t) starting from the Fock states with n0 = 0 . . . 7 respectively. Same color code as in Fig.2a. Dotted black lines are theoretical.

FIG. 4 :

 4 Photon number n

Acknowledgements This work was supported by the Agence Nationale pour la Recherche (ANR), by the Japan Science and Technology Agency (JST), and by the EU under the IP project SCALA. S.D. is funded by the Délégation Générale à l'Armement (DGA).