J.-P Thirion
email: jean-philippe.thirion@inria.fr

Image matching as a diffusion process: an analogy with Maxwell's demons

Keywords: deformable model, elastic matching, image sequence analysis, inter-patient registration, non-rigid matching

In this paper, we present the concept of diffusing models to perform image-to-image matching. Having two images to match, the main idea is to consider the objects boundaries in one image as semi-permeable membranes and to let the other image, considered as a deformable grid model, diffuse through these interfaces, by the action of effectors situated within the membranes. We illustrate this concept by an analogy with Maxwell's demons. We show that this concept relates to more traditional ones, based on attraction, with an intermediate step being optical flow techniques. We use the concept of diffusing models to derive three different non-rigid matching algorithms, one using all the intensity levels in the static image, one using only contour points, and a last one operating on already segmented images. Finally, we present results with synthesized deformations and real medical images, with applications to heart motion tracking and three-dimensional inter-patients matching.

INTRODUCTION

Many concepts of thermodynamics have been fruitfully applied in the field of information theory and, more specifically, to image processing. A recent example is the application of mutual entropy minimization techniques (see [START_REF] Viola | Alignment by maximization of mutual information[END_REF][START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF] to the matching of medical images acquired with different modalities.

Another instructive example is anisotropic filtering (see [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Catté | Image selective smoothing and edge detection by nonlinear diffusion[END_REF][START_REF] Gerig | Nonlinear anisotropic filtering of MRI data[END_REF][START_REF] Kimia | Geometric heat equation and nonlinear diffusion of shapes and images[END_REF], now commonly used as a preprocessing step for medical images: it can be shown that the application of a Gaussian filter of parameter σ is equivalent to the diffusion of heat in a homogeneous material for a time duration directly related to σ . Heat propagation is a new way of looking at Gaussian filtering: variations of this concept, for example by using the object boundaries in the image to create inhomogeneities with respect to heat propagation, have led to a large variety of new, non-linear algorithms for image filtering, even if convergence proofs have not always been established. Anisotropic filtering is not a new technique, but a new exciting way to consider a problem as ancient as image processing.

Extending the presentations in [START_REF] Thirion | Fast non-rigid matching of 3D medical images[END_REF][START_REF] Thirion | Non-rigid matching using demons[END_REF], we propose an original viewpoint for image-to-image matching, also based on an analogy with thermodynamic concepts (see Figure 1). Having two images to match, the main idea is to consider the object boundaries in one image as semipermeable membranes and to let the other image, considered as a deformable grid model, diffuse through these interfaces, by the action of effectors situated within the membranes. We illustrate this concept by an analogy with Maxwell's demons and we contrast it with more conventional viewpoints such as deformable models based on attraction. As we will see, diffusing models rely mainly on the notion of polarity (inside-outside), while attraction relies on the notion of distance, but mixed models can be devised, and we will also see how to translate the diffusing model concept into the concept of attraction, with an intermediate step corresponding to optical flow methods. First, we recall some existing techniques concerning nonrigid matching, such as 'snakes-based' methods, similarity maximization techniques and optical flow. Then we present and detail our concept of diffusing models, illustrated by its application to a simplified case, when two identical discs are matched rigidly (demon 0). We present a general iterative scheme to implement diffusing models, with several possible variants (demons 1, 2 and 3) which illustrate how the concept can be applied to generate new matching algorithms. Finally, we present experimental results with several medical applications such as the tracking of deformable organs or three-dimensional (3-D) inter-patients matching.

NON-RIGID MATCHING TECHNIQUES

Matching is an essential task for many computer vision applications. A clear definition can be given for rigid or articulated bodies: it is to recover rigid displacements of rigid parts. The task is much more complicated for deformable objects, with plastic or elastic deformations. In this case, there is not a single definition of an ideal optimal match, but as many definitions as practical applications. Each time, one has to define precisely the set of deformations T which are explored (rigid, affine, spline, free-form etc), and the type of features which are used (points, curves, surfaces, intensities etc).

However, classifying all non-rigid matching techniques with a single metric is made impossible in practice by the numerous works recently performed in that domain for a large number of different goals. We only describe a few of these methods here.

The concept of attraction

A widespread (and very intuitive) way to consider matching is based on an analogy with attraction. One example is gravitation: a point P of the deformable model M is attracted by all the points P in S which are similar. For example, let K (P, P ) be a similarity criterion, and D(P, P ) a function of the distance (not necessarily Euclidean), the force f on P, induced by the attraction of all the points of S can be f (P) = P ∈S K (P, P ) D(P, P ) P P .

(1)

M is deformed according to these forces, and also according to smoothness constraints internal to M. Such a method is computationally too expensive [O(n 2 )]. A less expensive method is to retain in the computation only one principal attractive point of S: a point P of M is attracted by the point P of S which is the 'closest and most similar' to P.

We perceive in this definition that a balance between being 'close' and being 'similar' must be determined. For example, in some methods, only contour points are used, and contour points in S are equally similar to contour points in M, making proximity the most important factor (see Figure 2). This is the case for iterative closest point (ICP) methods; see [START_REF] Besl | A method for registration of 3-d shapes[END_REF] and [START_REF] Zhang | Iterative Point Matching for Registration of Freeform Curves[END_REF] for the rigid case. This is also a basic assumption in the 'snakes' method introduced by [START_REF] Kass | Snakes: active contour models[END_REF]; see also [START_REF] Blake | Active Vision[END_REF] for a review of these techniques.

Distance and similarity can be married in an elegant way. For example in [START_REF] Feldmar | Locally affine registration of freeform surfaces[END_REF], the points have n attributes (or parameters), defined using differential geometry: these attributes are invariants associated with the points, such as principal curvatures. The coordinates x, y, z of the points are also considered as parameters, therefore each point of S is represented by a single point in an (n + 3)dimensional parametric space, where spatial coordinates and differential invariants are mixed. Then only distance is needed, like in ICP, but in the (n + 3)D parametric space: it can be the Euclidean distance, or the Mahalanobis distance if uncertainty is evaluated as well (see [START_REF] Ayache | Artificial Vision for Mobile Robots: Stereo Vision and Multisensory Perception[END_REF].

Similarity can also be introduced into snakes-based techniques. For example in [START_REF] Benayoun | Dense non-rigid motion estimation in sequences of 3D images using differential constraints[END_REF], a third term is added to the energy function to be minimized, corresponding to the similarity of differential geometry properties attached to those points. This helps to reduce the well known 'aperture problem' (see [START_REF] Faugeras | Three-Dimensional Computer Vision: a Geometric Viewpoint[END_REF], which states that it is easier to recover displacements normal to the contour rather than along it.

In contrast to ICP or snakes, correlation techniques place more importance on similarity than on distance. The point P of M is attracted by P of S which maximizes a correlation function K (P, P ). But a distance is also used, which is 1 if P is within a neighbourhood of T i (P) (where T i is the current estimated deformation) and +∞ if P is out of the correlation window.

As described, these methods are still expensive because at each iteration and for each P, a large number of possible attractive points P in S are to be considered. In practice, there are algorithmic solutions to reduce the complexity of finding this 'closest and most similar' point. It can be the use of KD trees, like in [START_REF] Besl | A method for registration of 3-d shapes[END_REF] or [START_REF] Feldmar | Locally affine registration of freeform surfaces[END_REF]. It can be by reducing the set of feature points, like using the surface of organs (see [START_REF] Davatzikos | Nonlinear registration of brain images using deformable models[END_REF][START_REF] Thompson | A surface-based technique for warping three-dimensional images of the brain[END_REF], or an even more compact representation such as crest lines (see [START_REF] Declerck | Automatic retrieval of anatomical structures in 3D medical images[END_REF] or individual points (see [START_REF] Thirion | Extremal points: definition and application to 3D image regist ration[END_REF][START_REF] Rohr | Point-based elastic registration of medical image data using approximating thin-plate splines[END_REF]. It can be the use of pre-computed distance maps, such as for chamfer matching techniques (see [START_REF] Borgefors | Hierarchical chamfer matching: a parametric edge matching algorithm[END_REF][START_REF] Malandain | Improving registration of 3-D medical images using a mechanical based method[END_REF].

An attractive point can also be considered to be a minimum of a potential field (an 'attractor'): by differentiating this potential field a local expression of the forces is obtained. This can be done for snakes [using, for example, the gradient, see [START_REF] Kichenassamy | Gradient flows and geometric active contour models[END_REF]], and also for correlation techniques. These differentiations give forces which are likely to be directed toward the 'closest most similar' point, but it is only a local approximation.

As we can see, the concept of attraction has inspired the majority of existing matching methods. Distance and similarity are central in attraction, while polarity is accessory and seldom used. We note however that there exist some works that make use of polarity-in [START_REF] Radeva | A snake for model-based segmentation[END_REF] polarity is used to discard the forces of contours whose normals are directed away from their closest point. In the work of [START_REF] Chakraborty | Deformable boundary finding in medical images by integrating gradient and region information[END_REF], regional information is used in the form of a third term to be minimized along with internal and external forces. This third term corresponds to the integral over the model interior of the points which are unlikely to be interior points. One can refer to the work of [START_REF] Ronfard | Region-based strategies for active contour models[END_REF] for a method mostly based on polarity: we will discuss this last method in detail later on.

Optical flow methods

A special kind of method is optical flow, which is used to find small deformations in temporal sequences of images (see [START_REF] Horn | Determining optical flow[END_REF][START_REF] Aggarwal | On the computation of motion from sequences of images-a review[END_REF][START_REF] Barron | System and experiment: performance of optical flow techniques[END_REF]. At a given point P, let s be the intensity function in S and m the intensity in M (see Figure 3). The basic hypothesis of optical flow is to consider that the intensity of a moving object is constant with time, which gives, for small displacements, the optical flow equation

v • ∇s = m -s. ( 2
)
This constraint is not sufficient to define the velocity v; see for example [START_REF] Simoncelli | Probability distributions of optical flow[END_REF]. One solution is to regularize the problem to obtain the local values of v. Another solution is to consider that the end point of v is the closest point of the hypersurface m, with respect to spatial (x, y, z) translations (see Figure 3), which leads to Equation ( 3

): v = (m -s) ∇s ( ∇s) 2 . ( 3
)
This equation is unstable for small values of ∇s, leading to infinite values for v. Ideally, the expression should be close to zero for small ∇s. A solution is to multiply Equation (3) by (( ∇s)

2 /(( ∇s) 2 + (ms) 2 )), which gives Equation (4):

v = (m -s) ∇s ( ∇s) 2 + (m -s) 2 (4) or, v = 0 if ( ∇s) 2 + (m -s) 2 < .
With this expression, the optical flow can be calculated in two steps: first compute the instantaneous optical flow for every point in S, then regularize the deformation field.

In optical flow, v is considered to be a velocity because the images are two successive time frames: v is the displacement during the time interval between the two image frames. In fact, when comparing images of two different patients, there is no such temporal consideration and it is more general to consider v as being simply a displacement.

DIFFUSING MODELS

We now present the concept of diffusing models, with a parallel with Maxwell's demons.

Maxwell's demons

The concept of demons was introduced in the 19th century by Maxwell to illustrate a paradox of thermodynamics (see Figure 4). Assume a gas composed of a mix of two types of particles a and b a , and separated by a semi-permeable membrane. Assume also, that this membrane contains a set of 'demons', which are able to distinguish between the two types of particles, and allow particles of type a only to diffuse to side A of the membrane and particles of type b only to diffuse to the other side B. At the end, A contains only particles a, and B particles b. This corresponds to a decrease of entropy, in contradiction with the second principle of thermodynamics. The paradox was solved because the demons generate a greater amount of entropy to recognize the particles; thus, the total entropy of the system has increased.

Demons for image processing

Let us see how to apply demons to image matching. We want to match a model image M with a scene image S: M must be deformed to resemble S as much as possible. We assume (see again Figure 1) that the contour of an object O in S is a membrane, and we scatter our demons along this contour. We assume also that we are able to determine locally, for each contour point in S, a vector perpendicular to this contour and oriented from the inside of the object to the outside (for example, the gradient of the image S can be used). We assume also that M is a deformable grid, whose vertices are particles which can be classified as 'inside' or 'outside' particles. The rigidity of the deformable grid M is determined by the relations between these particles, and different behaviours, from totally rigid to totally free form, can be obtained by varying these relations. We will see later on that many types of deformations can be used. We call M a diffusing model, and we give an informal definition for demons:

a Originally a and b were cool and hot particles. Polarity (inside-outside) is central for diffusing models, and distance is optional.

Diffusing models: a simplified example

To illustrate this principle, we consider the simple case of two images displaying the same disc, and we restrict ourselves to rigid transformations. We consider how this problem is tackled having in mind the concept of attraction as well as of a diffusing model (see Figure 5).

• For attraction, it is natural to sample regularly the model disc boundary (like in snakes-based techniques), and to apply forces directed toward the closest point of the scene circle, which contains the attractive points: distance is central. • For a diffusing model, we regularly sample the scene disc boundary, each point being a demon. The force of the demon is oriented from the inside to the outside of the disc if the corresponding model point is labelled 'outside', and the other way if the label is 'inside': polarity is central.

Both methods are iterative: at each iteration, the motion created by all the elementary forces is applied to the model. Again to simplify the comparison, the force magnitudes induced by the demons and the attractive points are constant and equal at each iteration, but decrease with each iteration, in order to allow convergence. In Figure 5, the top three images represent three iterations using attraction, and the bottom three images are three iterations using a diffusing model. We can see that the individual forces are clearly different between both methods, and that this is not due to a change of the reference system, because the transformation is rigid. Figure 6 presents the next three iterations for attractive points (top) and diffusion (bottom). The two methods converge toward a similar result (but this is not always the case) and for small deformations, the forces also become similar.

We note that if the two discs do not overlap initially, then the demon-based method does not apply, whereas attraction still works, but there are other examples where attraction models can get stuck into local minima (see Figure 7 for an example) while diffusing models can still give the expected solution.

It is remarkable to observe that a model solely based on polarity can lead to matching.

Related works

From all the methods that we have considered, the closest work in its spirit is probably the 'anticipating snakes' of [START_REF] Ronfard | Region-based strategies for active contour models[END_REF] used for two-dimensional (2-D) image segmentation (see Figure 8). The principle of the method is to consider a deformable contour where the external forces are replaced by 'region-based' forces, relying mainly on polarity, and which fit perfectly to our definition of demons. Anticipating snakes is distinct from diffusing model though, mainly because what is considered in [START_REF] Ronfard | Region-based strategies for active contour models[END_REF] is a deformable contour evolving in a static image instead of a deformable image grid diffusing through a static contour. One visible consequence is that forces are applied to a deformable contour in [START_REF] Ronfard | Region-based strategies for active contour models[END_REF], and to a deformable grid in our method. Along iterations, external forces are also always applied in the same location of the deformable contour in anticipating snakes while they are applied on successive locations of the deformable grid in diffusing models.

With respect to medical applications, the most related works are [START_REF] Bajcsy | Multiresolution elastic matching[END_REF], [START_REF] Gee | Elastically deforming 3D atlas to match anatomical brain images[END_REF] or Christensen et al. (1994b) where a global image, considered as a 2-D or 3-D grid, is also deforming. In these methods, the external forces exerted on the grid are traditional ones, based for example on the derivatives of a cross-correlation similarity measure.

It seems to us that the latter methods could be studied from the viewpoint of diffusing models, although it is still unclear to us how to do so. However, we are able to show that, to some extent, optical flow can be considered as an intermediate concept between diffusing models and attraction.

Optical flow behaves like diffusing models

Through each point P of a scene image S where ∇s = 0 is an iso-contour s = I , where I = s(P) is constant. This isocontour is the interface between the inside regions s < I and the outside regions s > I . Comparing the intensities of the model image M with I also gives an automatic way to label the points of M 'inside' or 'outside'.

The displacement v [Equation ( 4)] is comparable to the application of an elementary force f during one iteration step, whose direction is the same as ∇s and whose orientation b is according to -∇s if m < s, that is m < I and according to + ∇s when m > I . In other words, the force f pushes a point P of M toward the outside (that is, ∇s) when P is labelled outside (m > I ), and toward the inside (that is -∇s), when P is labelled inside (m < I ), which is exactly our definition of a demon.

Optical flow also behaves like deformable models

based on attraction Consider again Figure 3: v is the shortest spatial displacement that brings the point [P, m(P)] of the hypersurface corresponding to M, into the hyper-plane coming through [P, s(P)], with normal (-∇s, 1), which is the best local approximation of the hypersurface S. If only local information are available [that is m(P), s(P) and ∇s(P)], it is legitimate to assume that P = P + v is the point of S closest to P, and having the same intensity [m(P)]: P is the closest point in S having the same intensity than P in M, which is by definition the behaviour of attractive points.

IMPLEMENTATIONS DERIVED FROM THE CONCEPT OF DIFFUSING MODELS

We propose a general scheme with several possible variants, leading to different implementations of diffusing models.

Matching as an iterative process

Similarly to deformable models based on attraction, diffusing models require an iterative scheme. The aim is to find a final transform T ∈ T (where T is the set of allowed deformations), between the space M of the model image M and the space S of the scene image S. T is the final evolution of a family of transforms {T 0 , T 1 , . . . T i , . . .} ⊂ T . At each step, the deformed version T i (M) of the model M becomes T i+1 (M), constrained by 'internal' forces f int created by the relations between the model points, and 'external' forces f ext , created by the interactions between T i (M) and S. We can also change the reference system and process in M the interactions between M and T -1 i (S). Both of the simplified examples described previously (diffusing model and attraction) illustrate an iterative scheme: in these cases, T is the group of rigid transforms and T i+1 is δT i • T i where δT i is the residual motion created by the set of elementary demons or attraction forces.

b We distinguish between direction and orientation: a straight line has a direction. An oriented line, a half line or a vector has an orientation as well as a direction, which is the way it is headed.

Compute Ti+1 from Ti and {f}

Compute the forces {f} between Ti(M) and S Pre-compute the demons Ds In the case of diffusing models, the first step is the precomputation of the set of demons D s , extracted from S, and the second step is an iterative estimation of the deformation T , from the model space M to the scene space S (Figure 9).

Extracting the demons from S

The demons set D s is extracted from the scene image S. D s can be the whole image grid (one demon per pixel or voxel), such as in [START_REF] Thirion | Fast non-rigid matching of 3D medical images[END_REF]: in that case, the interfaces at each points are the iso-contours. D s can be also restricted to the contour points in S, such as in the simplified example (demons 0), or to points automatically extracted by edge detection methods. The information attached to each demon can be:

• its spatial position P in S (possibly sub-pixel); • a direction d from the inside to the outside [generally based on the gradient ∇s(P)]; • the current displacement from S to M at that point: d = P P , where P = T -1 i (P); • information about the interface, such as the intensity at that location s(P).

The iterative part

We start with an initial deformation T 0 (i.e. identity). At iteration i, we have a current estimated transform T i , and each iteration is composed of two steps:

(i) For each demon P ∈ D s , compute the associated elementary demon force f i (P) which depends on the demon direction d s at point P and on the polarity of M at point T -1 i (P). (ii) Compute T i+1 from T i and from all the elementary demons forces { f i (P), P ∈ D s }.

In some cases, step (ii) can be decomposed into two steps: first compute an elementary deformation δT i from the elementary forces f i , then T i+1 = δT i • T i .

Possible variants

Useful variants of the scheme above can be obtained by varying:

(i) the selection of the demon positions D s (whole image grid, contours points etc). (ii) the space of deformations T (rigid, affine, spline, free form etc). (iii) the interpolation method which gives the values of M for the non-integer positions T -1 i (P) (linear, spline, sinc etc). (iv) the formula giving the force f of a demon (constant magnitude, gradient-based, optical flow like etc).

Thus, the simplified diffusing model example (demons 0) uses: (i) sample points of the disc contour of S for D s ; (ii) rigid transforms for T ; (iii) no interpolation because m(T -1 i (P)) is analytically defined and (iv) constant magnitude forces. In the following, in addition to the simplified case, we examine three other examples of diffusing models, illustrating different variants of the general scheme.

Demons 1: a complete grid of demons

The following method is especially suitable for 3-D medical image analysis (see [START_REF] Thirion | Fast non-rigid matching of 3D medical images[END_REF]. The general scheme is modified as follows.

(i) All pixels of S where ∇s = 0 are selected to be demons (D s = S). (ii) T is a free form deformation, that is, for each demon P we store the current elementary displacement d(P). To get a regular displacement field, a Gaussian filter with a given σ is applied to the whole field at each iteration.

A theoretical study about how smoothing deformation fields can be used for regularization can be found in [START_REF] Anandan | A computational framework and an algorithm for the measurement of visual motion[END_REF]. (iii) m(P ), where P = P + d(P) is estimated using trilinear interpolation in M. (iv) the demon force is given by optical flow [ v of Equation (4)], using only the information provided with each point: m(P ), s(P) and ∇s(P). More precisely, we do not define the force, but the result of the application of a force during one iteration step, which is a displacement:

d = -v.
To make the algorithm faster (each voxel is considered at each iteration) and more robust with respect to initial positions, we adopt a multi-scale scheme, starting with a large number of iterations at coarse scales, down to a few iterations at the finest scale. More precisely, we use a pyramid approach were, at each scale, 1 8 of the voxels and four times the number of iterations are used than that in the immediately finer scale, which means that the computation time for the rest of the pyramid is equivalent to the time spent at the finest scale.

This gives us results qualitatively similar to other deformable grid methods such as those of [START_REF] Bajcsy | Multiresolution elastic matching[END_REF] or Christensen et al. (1994b). However, it is hard to assess quantitative differences, firstly because the definition of an 'ideal' inter-patients matching remains open to interpretation (see [START_REF] Thirion | Cross validation of three inter-patients matching methods[END_REF].

With respect to computation time, our method is particularly time efficient, when compared to the results described by [START_REF] Christensen | Individualizing neuro-anatomical atlases using a massively parallel computer[END_REF] with the use of a massively parallel computer. Besides the computational aspect of finding similar points, the method presented in Christensen (1994b[START_REF] Christensen | Individualizing neuro-anatomical atlases using a massively parallel computer[END_REF] is also computationally expensive because it models exactly the physical behaviour of a viscous fluid. [START_REF] Bro-Nielsen | Fast fluid registration of medical images[END_REF] have shown (and in 2-D) that those computations can be implemented using a convolution filter applied to the deformation field, leading to an improvement of one order of magnitude in speed. Another interesting result is that they demonstrate theoretically that the Gaussian filtering used in our method to regularize the deformation field approximates linear elasticity.

The algorithm demons 1 could have been derived from optical flow concepts too: to some extent, demons 1 can be considered as a particular kind of multi-scale, iterative version of optical flow. We derived it from diffusing models, which emphasize how optical flow can be treated as a diffusion process. This reveals possible changes in the algorithm, for example in the expression of the forces (see demons 2 and 3 below) which have less and less to do with the optical flow equation. We note however that using the optical flow formula for demons leads to particularly good results for the experiments presented in the last section.

Demons 2: demons in contours only

In matching methods, contour points are generally more important than other points in the image. Furthermore, using only the contour points of S leads to a faster algorithm. We used, for the possible variants:

(i) the Canny-Deriche edge detector (see [START_REF] Deriche | Using Canny's criteria to derive a recursively implemented optimal edge detector[END_REF][START_REF] Monga | 3D edge detection by recursive filtering and edge tracking[END_REF]) to extract D s from S; (ii) global transforms (rigid, affine) for T , determined by least squares [like in [START_REF] Besl | A method for registration of 3-d shapes[END_REF]].

It could be improved by extended Kalman filtering (EKF; see [START_REF] Ayache | Artificial Vision for Mobile Robots: Stereo Vision and Multisensory Perception[END_REF] to reject outliers. We have also tested more flexible deformations using warping techniques (see [START_REF] Bookstein | A feature space for derivatives of deformations[END_REF][START_REF] Szeliski | Matching 3-d anatomical surfaces with non-rigid volumetric deformations[END_REF][START_REF] Szeliski | Hierarchical spline-based image registration[END_REF][START_REF] Declerck | Automatic retrieval of anatomical structures in 3D medical images[END_REF] (iii) (tri-)linear interpolation;

(iv) the demon behaviour is detailed below.

We now describe the demon force computation in this case. As only a subset D s of S is considered, more sophisticated demons can be designed than in the previous method relying on the optical flow formula. Let P be a contour point in S, of intensity s(P), and n = ∇(P)/|| ∇(P)|| be the oriented normal of the contour (from inside to outside). By interpolation, we compute the value s out = s(P + k n) and s in = s(Pk n) (where k is a constant, possibly 1).

P is therefore the interface between a region whose estimated intensity is s in and a region whose intensity is s out . The demon force is

f (P) = K s in ,s out (m(P )) n ( 5
)
where P = T -1 i (P) and K s in ,s out (m) is a function that we now detail. By definition, the demon pushes inside (K < 0) when P is labelled 'inside' in M, and outside (K > 0) when P is 'outside' in M.

We have considered functions K where the probability of being 'inside' is very high if m(P ) ≤ s in (then we impose K = -k), and the probability of being 'outside' is very high for m(P ) ≥ s out (then K = k). We have used the simple piecewise linear formula of Figure 10 for K (m).

We have considered another type of function K , where when m is too different from s in and s out , P is neither considered 'inside' nor 'outside', but simply 'irrelevant' (K = 0). We assume that the intensities of inside points are s in , and outside points are s out , but that the measure s is corrupted with Gaussian noise. We use therefore for K the sum of two Gaussian distributions K = K -+ K + (see Figure 11); K -(m) is negative and centred at s in , with K -(s in ) = -k, and K + (m) is positive, centred at s out , with K + (s out ) = k.

We still have two free parameters (σ -and σ + ), which can be related to the level of noise in the images S and M. Anisotropic diffusion filters mentioned in the introduction could be used efficiently to pre-process S and M and reduce the noise level, while preserving the interfaces. Both of the K (m) functions described gave us good results and we kept the first one for the experiments in the present paper, but a more careful study is needed to determine when to use one or the other.

Demons 3: already segmented images

Here, the nature of the problem is different. We assume that the images S and M are already segmented, that is, for each point in image S or in image M, we have a label which refers to a given structure (for example a given organ in the case of medical images). Even though the images are segmented, we can still search for a geometric point to point correspondence between S and M. Within the general scheme, we use: (i) the set of interface points of S for D s , that is, points between adjacent voxels whose labels are different (Figure 12); (ii) deformations similar to method 2 (rigid, affine, warping); (iii) no interpolation, but the label of the closest voxel in M (using linear interpolation between labels is generally meaningless); (iv) demons with forces of constant magnitude, the orientations are defined according to the labels (see below).

We now detail the demons behaviour. Two adjacent voxels are linked by a segment of length 1, √ 2, √ 3 depending on their connectivity (6, 18, 26) (see Figure 12). Let A and B be two adjacent voxels whose labels s in = s(A) and s out = s(B) are different, and d = AB. We create a demon in P = (A + B)/2.

The 'force' f (P) of the demon is determined by s in , s out and the label m(P ), where For a given iteration, k is a constant magnitude, which can vary at each iteration depending on the convergence strategy. For example, k can be set to k 0 voxel size and decreased linearly to zero during iterations. Accordingly, the smoothing parameter σ can be decreased from σ 0 to 0, to ensure the convergence of the algorithm. The ratio k 0 /σ 0 depends on the application: it controls the rigidity of the diffusing model. The advantage of using pre-segmented images is that the interface points are labelled (s in , s out ), and interact with the structures of M which are either s in or s out , but do not interact at all with other structures.

P = T -1 i (P), then • if m = s in and m = s out , f = 0 • if m = s in , f = -k d • if m = s out , f = k d.
For automatic registration with a digital anatomical atlas, such a method is very suitable: a physician can spend a lot of time labelling all the voxels in the image S of a reference patient (see Figure 26), then a method using an hybrid of demons 2 and 3 can be applied to automatically match and label, at the same time, the image M of a new patient.

'Bijectivity' of the 3-D deformation

We do not enter into the debate of what bijectivity means in the case of discrete deformation fields, but one important constraint for 3-D deformations is the possibility of inverting these deformations. With the previously described implementations, there is no guarantee that this inversion can be performed [a problem which is raised in [START_REF] Bro-Nielsen | Fast fluid registration of medical images[END_REF]]. The bijectivity property is ensured to some extent in physically based techniques [such as those of Christensen et al. (1994a) or [START_REF] Bro-Nielsen | Fast fluid registration of medical images[END_REF]] by a local control of the Jacobian (the determinant |J | of the Jacobian matrix of the deformation), which is forced to be positive (|J | > 0). This type of bijectivity should not be confused with the 'one-to-one' correspondence of viscous fluids (see [START_REF] Christensen | Individualizing neuro-anatomical atlases using a massively parallel computer[END_REF], corresponding to incompressible matter (|J | = 1) which might be useful for specific types of problems.

However, the elastic fluid is probably a better model for inter-patients matching than viscous fluid a although there a Firstly, the two patients are likely to have different sizes. are no grounds to relate these two mechanical models to anatomical differences. In the following applications we will be mainly concerned the elastic deformations.

Maintaining a constraint such that |J | > 0 (elasticity) is also time consuming. To simply verify that a discrete 3-D deformation field has a positive Jacobian, one can consider each cubic cell (or 8-cell) composed of 8 voxels, and verify that the deformed 8-cell is still positively oriented. One way to do this is to decompose the 8-cell into five tetrahedra, and verify that the deformed tetrahedra are still positively oriented (by computing the determinant of the matrix formed by three edges), which means about 100 floating-point operations per voxel and per iteration, which would make this measurement a bottle neck of our method, without even solving the problem (which is to impose a positive Jacobian).

We have implemented a simpler solution (see also [START_REF] Burr | A dynamic model for image registration[END_REF]: at each iteration, we compute the direct deformation T 12 (from I 1 to I 2 , for a regular grid in I 1 ), and the reverse deformation T 21 (from I 2 to I 1 , for a regular grid in I 2 ), and also the residual deformation R = T 21 • T 12 . We then compensate equally for this residual deformation in T 12 and T 21 , by removing 'half' of R from each deformation. More precisely (Figure 15), a point P in image I 1 is projected into P = T 12 (P) in I 2 , and, by interpolating T 21 at the floating Hence the residual vectors P P constitute a vector field for a regular grid in I 1 . 1/2P P is directly subtracted from the vector field T 12 and the other half 1/2P P , is projected into I 2 by using T 21 and added to it, to give a new T 21 . After this operation, T 21 • T 12 is very close to identity, and we have verified experimentally that this residual is less than a voxel on average, when this operation is performed at each iteration of the diffusing model.

This method requires only twice the computation of the one way transformation plus only about 20 additional floating point operations per voxel per iteration. Besides bijectivity, the advantage of this method is that it provides at the end not only the direct transformation T 12 , but also an inverse transformation T 21 very similar to the 'ideal' T -1 12 . Experimentally, this bijective implementation has proved to be much more robust than the simple one for anatomically homologous organs (for example right hands with five fingers each). It tends to preserve the morphology of each specimen. For example, in the temporal lobes of the brain, the number of gyri for different patients can be different, but are still preserved: only low frequency morphometrical differences are compensated for, which allows for the comparison of local morphological differences (see [START_REF] Guimond | Automatic MRI database exploration and applications[END_REF].

Partial conclusion

To some extent, diffusing models (principally demons 1) can be connected back to more classical concepts of matching. However, the concept can serve to foster the exploration of new matching algorithms by the many possible options. Diffusing models are not constrained to use forces derived from a known potential field or from a similarity measure; the concept provides a global understanding of how things are working: the deformable model diffuses through the scene interfaces in order to have more points of the model with polarities or types similar to those at corresponding places in the scene.

EXPERIMENTS

We now present some experimental results of non-rigid matching obtained with implementations 1 and 2 of diffusing models, using synthesized and real data. We have used the same set of parameters for all the experiments: σ = 1 for the Gaussian filtering, four iterations at the finer scale of the pyramid, and four levels of resolution with respect to the multi-scale processing.

Experiments with simulated deformations

First, we present some results with a synthetic deformation in 2-D. Figure 13 shows a typical slice from a magnetic resonance (MR) brain image, and the same image after a deformation based on sinusoidal functions with a spatial period of 32 voxels.

We applied demons 1 to get the image of Figure 14 (CPU time, 30 s), where we have also superimposed the contours of the original image to verify visually the quality of the match, which is quite good, when compared with the large amount of deformation.

We have also successfully tested demons 2 (with warping for the deformations) on this image (see Figure 16). The quality of the match is also good, but slightly less satisfying than with demons 1. The difference, however, is very small.

The computation time is only 9 s in 2-D (three times faster, edge extraction included). This is a general result of our experiments: demons 2 is faster but less accurate than demons 1, because the density of demons is reduced in demons 2. The difference in quality is generally on the order of the voxel size. The choice of the method depends on a compromise between quality and computation time: none is better per se.

Lastly, we performed a test with the implementation demons 3, using warping for deformations. In this experiment only already labelled images are taken into account (the grey-level images are ignored, see Figure 17). This demonstrates that polarity (or labelling) is sufficient to perform successfully the non-rigid matching of two images.

3-D image sequences displaying cardiac motions

The two major drawbacks of non-rigid matching are reduced for image sequence analysis: the acquisition device being the same, the intensities of the scanned object do not vary much from frame to frame, and the position of the object in one frame is a good initialization for the next one.

Many deformable surface techniques have been proposed to study cardiac sequences, such as in [START_REF] Bardinet | Fitting of iso-surfaces using superquadrics and free-form deformations[END_REF], where a parametric surface is fit to the data and the deformation is analysed using modal decomposition. A very interesting study about the relevance of surface deformable models to heart motion analysis can be found in [START_REF] Shi | A unified framework to assess myocardial function from 4d images[END_REF], based on clinical experiments with dogs. However, these techniques require a segmentation of the object surface and the problem of extending the deformation to the whole volume remains to be solved.

In the following experiment, we present an example of the analysis of a CT scan of the heartbeat of a dog, acquired with the Mayo Clinic DSR (courtesy of Dr Richard Robb). To consider extreme conditions, we match directly the diastolic and systolic 3-D images. In Figure 18, we present two corresponding slices before elastic matching, and in Figure 19, the same slices after 3-D matching and resampling. With 3-D tools, we verified visually the quality of the match in all the other parts of the volumes. To get an idea of the smoothness of the non-rigid transform, we artificially added a regular grid tag to the image before deformation (see Figure 20). The computation time, for the matching of the two 100 3 voxels images was about 5 min on a DEC alpha workstation.

Figures 21 and22 present similar experiments for nuclear medicine 64 3 SPECT images (a human heart), with a CPU time of about 1 min. Figure 23 is a 3-D visualization of the same results.

We have also analysed whole cardiac gated sequences (eight 3-D images for a whole cardiac cycle). According to the physics of the SPECT acquisition, we have added two constraints: the cyclicity of the sequence, which is imposed by another set of iteration taking the whole sequence into account and performing temporal filtering, and the conservation of the radioactive matter. This last constraint states that, when the myocardium contracts, the density of radioactive material is increased, and thus the image intensity is correspondingly enhanced. This can be easily modelled within the design of the demons, by multiplying m(P i ) by the Jacobian of the transform, which represents the local variation of the density. The analysis of one whole time series takes only about 10 min. To conclude, we have an automatic, fast algorithm to estimate a deformation field between frames in 3-D image sequences, that does not necessitate a segmentation of the object, which is advantageous for SPECT images. A good point is that the physics of nuclear medicine acquisition can be taken into account in our scheme. The challenge now is to extract from this 4-D flow field the parameters that precisely characterize a pathology. Of course, these examples are given only to illustrate potential applications of diffusing models, but much harder work remains to be performed to validated the results from a medical viewpoint and to lead to useful clinical applications.

Inter-patients matching: the 3-D case

Matching the images of two different patients is very important for medical image applications. We present in Figure 24 slices from two different patients, extracted from their 3-D 256 × 256 × 128 MR images. We applied demons 1 and demons 2, but in 3-D, to the two blocks of data. The computation time, with the bijective computation, is only about 30 min on a Dec Alpha workstation, which is reasonably fast, with respect to the size of the data under consideration. Figure 25 presents the result of the matching: differences between the deformed patient and the reference patient is now hard to perceive. Only anatomical details which are topologically different between the two patients are still different. Again, demons 1 performs slightly better than demons 2, but at the expense of more computation time.

Because we have a 3-D completely labelled version of the reference image (courtesy of Dr Kikinis, Brigham and Women's Hospital, Boston, see Figure 27), we are able to apply the computed deformation field to it, and superimpose result on the image of the new patient (Figure 26), providing an automatic labelling of the patient. We note, however, that the inferred labels are mostly correct only where there are no local morphological differences between the two subjects. We hope to improve the automatic labelling with a more local, final pass of segmentation, taking the tissue types into account, such as in the segmentation presented by Wells III et al. (1995). Dr Benoît Dawant, at the University of Vanderbilt, TN, has used our method to compare the automatic labelling of images obtained by inference thanks to the deformation field that we can produce with several independent manual segmentation (see [START_REF] Dawant | Automatic 3d segmentation of internal structures of the head in MR images using a combination of similarity and free form transformations[END_REF].

A database of brain images

Finally, we present the results of the automatic non-rigid matching performed on a database of brain images. In order to perform this experiment, we have modified our method to estimate automatically a global bias and gain between the two images to be matched. We do this simply by incorporating into the iterative scheme an estimation of the two parameters based on the fitting of a line (with outliers rejection) to the transfer function between the two images, which gives, for each intensity in one image, the average intensity in the other image, based on the correspondence T i estimated at iteration i. We can see in Figure 28 a zoom on the central part of the brain in nine different patients, which shows both the differences in intensities and in shapes. Figure 29 shows the result of the automatic matching: after re-sampling, morphometrical differences between patients are hard to perceive, while local morphological differences still exist.

Figure 30 shows several slices of the reference image, and the average of the 10 images before and after matching with this reference image. The whole process is entirely automatic, and it takes about 5 min CPU time to match two 128 3 3-D images. It shows the quality of the matching, and also that the central parts of the brains are indeed more stable anatomically than the peripheral parts (which appear fuzzier).

Other medical applications

Our non-rigid matching method has been embedded into a more general system to automatically explore large databases of 3-D medical images in order to find similarities between patients, and to compute average representative specimens (see [START_REF] Guimond | Automatic MRI database exploration and applications[END_REF]. It is also being used for an ongoing large clinical study about Schizophrenia (E.C. project Biomorph) to measure and compare brain asymmetry between individuals. It has also been used for the computation of average activation maps in SPECT images (see [START_REF] Migneco | Spect memory activation studies thanks to non-rigid automated 3D image registration[END_REF]: similar ongoing studies are under way on a larger scale for fMRI and for PET activation images using competing non-rigid matching techniques (see [START_REF] Collins | Automatic 3D intersubject registration of MR volumetric data in standarized Talairach space[END_REF]. Finally, we have used our method to study and quantify the evolution of multiple sclerosis plaques in time series of 3-D brain MRI (see [START_REF] Thirion | Measuring lesion growth from 3D medical images[END_REF].

Validation with real images

The major problem of inter-patients matching is how to validate the different methods. In [START_REF] Thirion | Cross validation of three inter-patients matching methods[END_REF], we carefully compared three different inter-patients matching techniques with a database of CT scans of dry skulls: one matching method was based on the feature points defined by an anatomist, one was based on deformable 3-D crest lines, and the last technique was the demons-based method. The main results of this study were:

• the three matching methods gave mutually consistent results, with an average point to point distance of about 3 voxels; • this average distance was reduced to about 2 voxels when used to compute average patients' models; • the average influence of the choice of a reference patient on the resulting average model computations was less than 1 voxel. As mentioned, a study being performed at the University of Vanderbilt with our tool provides another independent validation of the application to automatic segmentation.

VIDEO DEMONSTRATIONS

Two video sequences illustrating the 3-D deformation field measured between two patients can be obtained from: http://www.inria.fr/epidaure/Gallery/atlas movie.html and atlas movie2.html.

CONCLUSION

We have proposed the concept of a diffusing model for image-to-image matching, and the related concept of demons to emphasize the role of polarity information in matching methods. A diffusing model corresponds to the 'diffusion' of a deformable model grid into the image of a scene, where the boundaries of the scanned objects are considered to be semi-permeable membranes, filtering the points of the model according to their polarity (inside or outside). We have also shown that optical flow can be considered as an intermediate step between diffusing models and more classical matching methods based on attraction.

Diffusing models and demons can be tools for examining matching methods from an original viewpoint, to generate ideas leading to new methods, or to improve existing ones by taking into account polarity, whose importance is generally underestimated. Considering matching as a diffusion process might also be a way to introduce new ideas from thermodynamics into image matching.

With respect to applications, we have shown different ways to apply these ideas to the medical field, such as organ tracking or 3-D inter-patients matching. Developing and testing such 3-D non-rigid matching tools is now a very important challenge for diagnosis, for surgery planning and control, and more generally, for a better understanding of the human anatomy and its variability.

Figure 1 .

 1 Figure1. Diffusing models: a deformed image, considered as a deformable grid, is diffusing through the contours of the objects in the static image, by the action of effectors, called demons, situated in these interfaces.

Figure 2 .

 2 Figure 2. Deformable model with attraction.

Figure 3 .

 3 Figure 3. Instantaneous velocity from image M to image S.

Figure 4 .

 4 Figure 4. Maxwell's demons and a mixed gas.

Figure 5 .

 5 Figure 5. Three iterations of a model based on attraction (top row) and a rigid diffusing model (bottom row). These examples are produced by actual implementations.

Figure 6 .

 6 Figure 6. Next three iterations (top row, attractive points; bottom row, demons).

Figure 7 .

 7 Figure 7. Example of problematic initializations: left, when the two objects to be matched do not overlap, diffusing models are inefficient. Right, with an attraction model that does not take polarity into account, and with forces decreasing with the distance, the model can get trapped in a local minimum.

Figure 8 .

 8 Figure 8. Two different uses of the concept of demons. Left, in 'anticipating snake' for image segmentation; right, in a diffusing model, used to perform image-to-image matching.

Figure 9 .

 9 Figure 9. Iterative scheme in the case of diffusing models.

Figure 10 .

 10 Figure 10. A first type of K functions.

Figure 11 .

 11 Figure 11. A second type of K functions.

Figure 12 .

 12 Figure 12. Demons in a labelled image.

Figure 13 .

 13 Figure 13. MR head slice: original (left), deformed (right).

Figure 14 .

 14 Figure 14. Deformations corrected with demons 1 (left), with the contours of the original image superimposed (right).

Figure 15 .

 15 Figure 15. The bijectivity is maintained by the computation of R = T 21 • T 12 and the redistribution of the residual at each iteration of the diffusing model.

Figure 16 .

 16 Figure 16. Original (left), corrected with demons 2 (right).

Figure 17 .

 17 Figure 17. Top left, original label image (I 1 ); top right, deformed label image (I 2 ). Bottom left, I 2 deformed toward I 1 using the implementation 'demons 3'; bottom right, deformed I 2 with a superimposition of I 1 contours.

Figure 18 .

 18 Figure 18. Corresponding diastolic and systolic slice before matching (dog).

Figure 19 .

 19 Figure 19. Corresponding diastolic and systolic slices after 3-D matching and re-sampling.

Figure 20 .

 20 Figure 20. Artificially tagged image (2-D grid of planes) to show the 3-D deformation.

Figure 21 .

 21 Figure 21. Corresponding diastolic and systolic slice before matching (SPECT, human).

Figure 22 .

 22 Figure 22. Corresponding diastolic and systolic slice after 3-D matching and re-sampling (SPECT, human).

Figure 23 .

 23 Figure 23. 3-D visualization of the left ventricle of the heart in SPECT images. Left is the diastole, middle is the systole, right is the surface of the diastole deformed toward the systole by the 3-D vector field obtained with demons 1.

Figure 24 .

 24 Figure 24. Two slices of two different patients (256 × 256 × 128 voxels).

Figure 25 .

 25 Figure 25. The two different patients after automatic matching.

Figure 26 .

 26 Figure 26. The deformation field is used to adapt the 3-D anatomical atlas presented in Figure 27 to the new patient. Top row: one slice of the MR image of the originally labelled subject along with the corresponding labelling. Bottom row: left, a slice of a new subject and right, the inferred labelling obtained by the application of the deformation field to the labelled reference subject. The inference is visually satisfying, up to the local morphological differences.

Figure 27 .

 27 Figure 27. A 3-D display of the segmented image of the brain, courtesy of Dr Ron Kikinis, and rendered by Gerard Subsol.

Figure 28 .

 28 Figure 28. One slice (out of a 128) of the original images of nine different patients. Shapes and intensities are very different.

Figure 29 .

 29 Figure29. The same nine patients after non-rigid matching, resampling and intensity correction. The computation is performed entirely in 3-D. Note that the morphometrical differences are compensated for, but not the local morphological differences. Our technique has been already applied, with qualitatively identical results, to more than a hundred cases.

Figure 30 .

 30 Figure 30. Left, five consecutive slices of the reference image; right, the average of 10 patients without matching; centre is the average of the 10 patients after matching and re-sampling. The fuzziness of the edges of the brain in the middle images express the remaining local morphological differences between the subjects, which have been reduced.

ACKNOWLEDGEMENTS

I wish to thank Nicholas Ayache, Michael Brady, Morten Bro Nielsen, Jérôme Declerck, Gérard Subsol, Xavier Pennec, Hervé Delingette, Sylvain Prima, Alexandre Guimond and the Epidaure people in general for stimulating discussions and specific comments about the ideas presented in this paper. Special thanks for Gregoire Malandain who implemented the Canny-Deriche-Monga edge detection, and to Allen Sanderson for his careful proof-reading of early versions of this paper.


Johan Montagnat
email: johan.montagnat@sophia.inria.fr

Herv Delingette
email: herve.delingette@sophia.inria.fr

Globally Constrained Deformable Models for 3D Object Reconstruction

Keywords:

To achieve geometric reconstruction from 3D datasets two complementary approaches have been widely used. On one hand the deformable model framework locally applies forces to t the data. On the other hand, the non-rigid registration framework computes a global transformation minimizing the distance between a template and the data. We rst show that applying a global transformation on a surface template, is equivalent to applying certain global forces on a deformable model. Second we propose a scheme which combines the registration and free-form deformation. This globally constrained deformation model allows us to control the amount of deformation from the reference shape with a single parameter. Finally, w e propose a general algorithm for performing model-based reconstruction in a robust and accurate manner. Examples on both range data and medical images are used to illustrate and validate the globally constrained deformation framework.

Introduction

In this paper we are interested in model-based reconstruction of data and the segmentation of images. Model-based approaches proved to be powerful for many reconstruction tasks by providing an a priori shape knowledge of the data to recover. We intend to build a robust reconstruction algorithm for the modeling of 2D and 3D data using respectively deformable contours, and deformable surface meshes.

The model-based reconstruction problem can be summarized in the following manner: let M be a contour or surface model of a given object consisting of a set of vertices fV i 2 IR d g, ( i = 1 : : : n ). Let D be a 3D dataset. The problem Preprint submitted to Elsevier Preprint 11 August 1998 is to nd a geometric transformation T such that T M] is an appropriate representation of the object in D. Two main model-based approaches can be found in the literature for the computation of T: free-form deformation and registration.

Free-form Deformation

With free-form deformation (FFD), the transformation T is only de ned at the model vertices. There is no a priori restriction on the transformation which c a n be seen as a set of independent displacement v ectors: 8V i 2 M ;! T V i ] 2 IR d .

The active contours originally proposed by Kass et al. 9] derive the transformation T from external forces representing data attachment and internal, or regularizing, forces. The external forces are devised to deform the model towards the data. The internal forces are based on the intrinsic geometric properties of the model used. The FFD scheme is illustrated in gure 1. Each model vertex V i is thus submitted to a weighted sum of internal, f int i , and external, f ext i , forces. We use a Newtonian law of motion to update each vertex, V i , position: m d 2 V i dt 2 = ; dV i dt + f int i + f ext i , where and are weights and is a damping coe cient. By discretizing time and space with nite di erences we get an explicit scheme:

V t+1 i = V t i + ( 1 ; )(V t i ; V t;1 i ) + f int i + f ext i (1)
where V t i is the location of vertex V i at time t. In this equation, the time step t does not appear. It is part of the coe cients , a n d . FFD is very exible since the only regularizing constraints are internal forces computed locally at each v ertex. It is also quite sensitive to noise and outliers. Section 4 shows FFD a example and associated problems.

The Registration Framework

Registration of a template consists of nding a global transformation T minimizing the distance between the transformed template and the dataset. T is constrained to belong to a transformation space T reg , usually a group for the composition operator: T 2 T reg 8P 2 IR d ;! T P] 2 IR d .

Typical registration methods iteratively compute the best transformation through least squares estimation until some convergence criterion is met. Besl 1] and Zhang 22] proposed an Iterative Closest Point algorithm (or ICP). The general registration scheme is illustrated in gure 2. As will be shown in section 5, registration provides a robust framework for shape recovery by restricting the set of allowed transformations. Yet, it su ers from limited shape variability o r computation instability problems.

Known Problems and Previous Work

For many shape-based reconstruction problems, the registration approach d o e s not have enough degrees of freedom (DOF) for the recovery of complex and high variable shapes, contrary to FFD which m a y not enforce the regularity o f the deformation to make the model robust to noise and data outliers. Internal forces, even computed at a large scale over the surface, are not equivalent to restricting the numberof DOF in a m o d e l .

Several approaches have been proposed in computer vision for controlling the amount of allowed deformation of a model. A second approach consists of including global parameters in the regularizing functionals. In 18], translation and rotation parameters are included in the internal energy to take, explicitly, i n to account rigid motions. In 5], Delingette introduces a rigidity parameter that controls the extent of the regularizing forces.

Globally constrained deformation

In this section we provide a link between computing a global transformation and applying an external force eld on a model M. W e then propose a uni ed framework encompassing both FFD and registration.

Registration as a Displacement Field

Let V t i bethe vertex position of M at time t. In the registration framework, we rst nd the besttransformation T t and nally apply the transformation at each vertex: V t+1 i = T t V t i ]. The displacement eld applied between two iterations is therefore T t ; I where I is the identity transformation: D = T t ; I = fT t V t i ] ; V t i g. To illustrate we 3. apply some external forces on a cross-shaped surface model and we plot the displacement vector at each model vertex. Figure 3(a) shows the displacement eld resulting from a rigid transformation while gure 3(c) shows the displacement eld of FFD.

In equation 1, we have seen that the displacement eld applied in the FFD is equal to:

fV t+1 i ; V t i g = f(1 ; )(V t i ; V t;1 i ) + f int i + f ext i g.
We can therefore consider that the registration method based on the ICP algorithm is equivalent to having a deformable model without inertial and internal forces ( = 1 , = 0 ) and submitted to the global force f global i = T t V t i ] ; V t i . The global force is an external force which leads to a deformation eld with fewer DOF than the usual external forces since it is constrained by the nature of the transformation T.

Globally Constrained Deformation

We have demonstrated the equivalence between registration-based deformation and the application of a global force to a deformable model. We now propose a constrained deformation scheme, where a deformable model is submitted to global, external, and internal forces. The purpose of this scheme is to have a computer-e cient deformable model with an easy control on the number of DOF. Our approach is to weight with a single parameter the in uence of the global forces versus the local forces (see gure 4):

V t+1 i = V t i + ( 1 ; )(V t i ; V t;1 i ) + f int i + f ext i + ( 1 ; )f global i (2)
We call , the locality parameter. It controls the number of DOF of the de-

Implementation

To motivate the introduction of the GCD framework, we perform experiments using our deformable model implementation based on simplex meshes. T h i s i s a discrete representation of models whose formalism applies in 2D (deformable contours or snakes) as well as 3D (deformable surface meshes).

Simplex meshes 6] are meshes with constant vertex adjacency. They allow a local description (at each vertex) of the shape they represent. They may model objects with almost any topology without dealing with surface parameterization problems. 2-simplex meshes are a natural extension of snakes in 3D and they provide a powerful framework to express regularizing constraints.

The simplex mesh framework is computationally very e cient since local force computation does not require a minimization step.

For internal forces, we are using a shape constraint that enforces local shape around a vertex up to a similarity. It provides an interesting behavior to models in parts where no dataset is available: they tend to keep locally the reference shape. Registration and FFD require the model vertices to match data points. To enforce robustness against outliers, we search for a data point to match vertex V i only within a limited distance of V i . We compute the external force f ext i as a vector directed along the normal direction proportional to the distance from V i to the dataset. [START_REF] Cootes | Active shape models, their training and application[END_REF] Free-Form Deformation Example FFD allows the model to deform with few restrictions. This enables the reconstruction of complex shapes such a s t h e v entricle shown in gure 5(d) starting from a very di erent template (here a 800 vertices sphere). However, in case of noisy or sparse data, FFD leads to an under-constrained system and the model is likely to lock o n u n wanted local energy minima. An under-constrained model is also likely to intersect itself.

Registration Examples

We implemented rigid transformations (6 DOF), similarities (7 DOF), and a ne transformations (12 DOF) using the closed-form expression proposed in 14]. We also investigated a cubic B-spline transformation evaluation (many DOF depending on the numberof knots) based on a gradient descent 4]. It appears that increasing the number of DOF improves accuracy of the t. We measure the goodness of the t by computing the accumulated distance between each vertex and its closest data point (see gure 7). In the case of cubic B-spline deformations, the number of DOF is su cient to get a close approximation of the data. However, the entailed deformations do not preserve the a priori geometry of the model (for instance, due to the initial position of the model, nose vertices migrated toward the left jaw).

In uence of the Degrees of Freedom

Computation Cost Issues

As can be seen in gure 7, increasing the numberof DOF causes an augmentation of the computation time. In particular, the B-spline transformation is far more expensive than other transformation classes. and a ne transformations. For the B-spline transformation, the LSE algorithm uses a gradient descent method (see 4] for details) which i s m uch more costly than the matching algorithm.

Numerical Stability Issues

Another issue is the stability of the method used. The number of equations needed is proportional to the number of DOF of the transformation. In the case of an a ne transformation, only 12 parameters are evaluated. With more DOF, a su cient number of equations, i.e. a minimum number of matches between model vertices and data points, must be found to stabilize the equation system. We encountered numerical stability problems with B-splines when data is noisy and no enough matches can befound.

A General Algorithm for Model-Based Segmentation

To o vercome problems encountered with the FFD and the registration framework we proposed a GCD scheme that provides an intermediate behavior between these two ends. In fact, we integrate the GCD in a "coarse to ne" algorithm that gradually increases the amount of deformation allowed during the deformation process.

Hierarchical Registration

In gure 6(e) we show a face reconstruction example using B-spline. The deformed surface is close to data points but the high numberof B-spline DOF combined with a rough initialization lead to the destruction of the geometrical correspondences during the deformation process. To achieve an accurate reconstruction with proper geometrical correspondences, we propose to gradually increase the number of DOF. By starting the transformation with few DOF, we can enforce the robustness of the reconstruction. Figure 8 shows the face registration result, when rst using rigid transformations, then similarity, then a ne, and nally B-spline transformations. We propose two strategies to set the low threshold. An absolute threshold value can be provided by the user as a percentage of the initial mesh size. A relative threshold is computed as a percentage of the initial displacement. After evolving a few iterations (three for instance), we compute the mean displacement and we set the threshold as a fraction of the obtained distance.

We made the face registration experiment with automatic thresholding: each deformation stage ends when the model displacement is lower than 0.1% of initial displacement. Figure 10(a,b) shows the resulting mesh. Figure 10(c) compares the manual and automatic threshold strategies. This reveals a signi cant diminution in the total computation time. More convergence results are given in section 7.2.3.

From Rigid Registration to Free-form Deformation

A natural idea is to extend the previous incremental process using GCD. To optimize computational cost and numerical stability we use transformations with few DOF such as rigid, similarity, or a ne transformations. The locality parameter allows us to tune the numberof DOF. GCD extends the registration framework in a computationally e cient manner since the search for the closest points has to be performed only once for the computation of the global and external forces. The resulting deformation scheme is illustrated in gure 11. We rst increase the number of DOF by making the global transformation evolve. We then gradually increase using the same displacement criterion. We can therefore have an automatic model-based reconstruction algorithm.

It is di cult to estimate the maximal value that the locality parameter should have (do we go up to FFD or do we stop for a maximal value?). With high values the model is more likely to t the data but it is also more sensitive t o noise. Thus the maximal value choice is a trade-o between reconstruction accuracy and data quality. We let the user decide.

λ = 0 Free-form Deformations 0 < λ < 1 λ = 1 λ = 0 λ = 0

Results

In this section, we present reconstruction and segmentation examples to validate our algorithm.

Qualitative Results

We consider a multi-modality face reconstruction example by tting a face template on an MR image of the head. Figure 12 Image segmentation is a key issue in many computer vision applications. We demonstrate here the ability of the GCD to automatically segment anatomical structures. For example, the liver is a soft organ that exhibits drastic interpatient v ariations. The model has to be deformable enough to take i n to account the shape variations. Abdominal organs have close grey-level values and they are di cult to isolate. Therefore the model has to beconstrained enough to berobust against noise and weak contours while deforming.

Figure 13 shows a liver segmentation example from an abdominal CT-scan. One slice of the CT-image can be seen in (a), the template (b) is rst registered then iteratively deforms by gradually increasing the locality parameter up to 35% (c). The result can be seen as the model trace in one slice of the data (d). Another application is the segmentation of blood vessels and aneurisms from angiographic images. For the segmentation of blood vessels, we use cylindrical models with a speci c axial constraint. The axial constraint is a class of global transformations that makes the surface model bendalong its axis. Details on the axial constraint computation can be found in 12]. Figure 14

Quantitative Results

The quantitative results given here validate the GCD accuracy for reconstruction and segmentation tasks.

Reconstruction accuracy

We ran the face reconstruction example proposed in section 5 using the GCD framework. The rigid, similar, and a ne registration stages are identical. For high DOF tting, we use GCD with a ne constraint instead of B-spline registration. Figure 15 We then t the mesh M ? on the original range data for a xed number of iterations to get model M 0 . Running the deformation process with di erent deformation schemes, the quality of registration is evaluated as the pointwise distance between the reconstructed model M 0 and M: d = P i kV i ; V 0 i k 2 .

A small value of d implies that vertices positions of M 0 are close to those of M and therefore the transformation is robust and accurate. Figure 16(c) shows the values of d for a rigid registration ( rst point), an a ne registration (second point), or a GCD (with a ne constraint) following a rst rigid t (other points). We get the best results at the intermediate values of the locality parameter whereas FFD or global (a ne) transformation lead to maximum distance values. This result can be interpreted in the following manner : if the model is too constrained, it cannot deform enough to t the dataset. On the contrary, i f i t h a s too many DOF, the surface vertices can be subject to large displacements during the deformation.

Convergence results

We used our algorithm to perform segmentation of the liver in images of an abdominal CT-scan database (an example of a particular image is given in gure 13). Figure 17 reproduces convergence results such as the one given in section 5 but at a larger scale. Figure 17(a) shows the number of iterations set manually (bars) for each deformation step. These thresholds have b e e n c hosen so that the deformation process properly segment any of our 30 liver images. They had to be high enough to allow convergence in any case. The solid line of gure 17(a) shows the mean number of iterations used by automatic thresholding and the minimal and maximal values.

Conclusion

We have introduced a general reconstruction framework that encompasses both deformable models and registration approaches. GCD is an e cient and simple algorithm for controlling the amount of deformation. It leads to good conservation of shape properties (such as curvature) during the deformation.

We h a ve shown that most accurate results are obtained by rst globally registering the template. As the model converges, we increase the number of DOF by using the GCD scheme. The evolution of the parameter allows us to gradually increase the amount of deformation from a ne transformation to free-form deformation.

Using deformable models we bene tfrom an a priori shape and geometrical information relevant throughout the deformation process. It is possible to obtain drastic shape variations of templates. By merging registration with the free-form deformation framework, we make the deformation process robust to noise and outliers while remaining computationally e cient. However, the total amount of deformation allowed is left to the user. By constraining the model, we reduce problems that can arise with FFD such as model self intersection but we do not strictly avoid them. The whole deformation algorithm is also very dependent on the closest point computation.

In the future, we plan to incorporate some additional statistical information in the reference model, in order to introduce meaningful deformation constraints. The statistical study of medical databases should allow for the introduction of statistically relevant deformation modes.

Fig. 1 .

 1 Fig. 1. The iterative s c heme for free-form deformation.

Fig. 2 .

 2 Fig. 2. The iterative s c heme for template registration.

Fig. 3 .

 3 Fig. 3. Di erent displacement elds on a cross-shaped surface.

Fig. 4 .

 4 Fig. 4. The iterative s c heme for globally constrained deformations.

Figure 5

 5 Figure 5 shows an example of FFD used to segment a left ventricle from an MR image.Figure 5(a) shows three orthogonal slices of the volumetric image (b) is the gradient image for the corresponding slices (c) shows the model inside the volumetric image initialized as a sphere it deforms toward image gradient points (d) is the resulting model after a re nement stage.

  Figure 5 shows an example of FFD used to segment a left ventricle from an MR image.Figure 5(a) shows three orthogonal slices of the volumetric image (b) is the gradient image for the corresponding slices (c) shows the model inside the volumetric image initialized as a sphere it deforms toward image gradient points (d) is the resulting model after a re nement stage.

Fig. 5 .

 5 Fig. 5. Free-form deformation example.

Figure 6

 6 Figure6shows a face reconstruction example using the four transformation classes presented above. The dataset is a cloud of 3D points acquired with a Cyberware scanner. The template (rendered surface) represents another face. The role of the template is to provide an a priori information on the shape to recover. Thus it is geometrically adapted (vertices are concentrated in high curvature areas).

Fig. 6 .

 6 Fig. 6. Face registration with di erent transformations.

Figure 7

 7 Figure 7 also compares the computation time of the four transformations used.In general, increasing the number of DOF increases the computation time since it involves the convergence of a minimization stage in the ICP algorithm. However, a ne transformations are slightly more e cient to compute than rigid ones (see 14] for details).

Fig. 7 .

 7 Fig. 7. Comparison between transformation classes.

Fig. 8 .

 8 Fig. 8. Registration with increasing numberof DOF.The remaining problem is to determine when to change the transformation class and therefore increase the number of DOF. Consider gure 9 which represents for each registration stage the total displacement of all model vertices between two successive iterations of the ICP algorithm. Due to the fast convergence of the ICP algorithm (notice the log-scale on the Y axis), the displacement is sharply decreasing. When the displacement is low enough the model does not evolve signi cantly anymore. Therefore we can set a low threshold to stop a deformation stage and increase the number of DOF by changing the transformation. This leads to good results with preservation of the model geometry as can beseen in gure 8(e).

Fig. 9 .

 9 Fig. 9. Displacement of model.

Fig. 10 .

 10 Fig. 10. Registration with automatic thresholding.

Fig. 11 .

 11 Fig. 11. General framework.

  (a) shows a face template. The model is embedded in the volume image (b) then registered and deformed toward the edge point using the GCD. The resulting surface is shown in (c). A control of the result (the model intersection with three orthogonal image slices) can beseen in (d).

Fig. 12 .

 12 Fig. 12. Face reconstruction from an MR image.

Fig. 13 .

 13 Fig. 13. Liver segmentation from contrast enhanced spiral CT-scan image.

  (a) shows three orthogonal slices of an angiography. Models are initialized around the aneurism and the vessels (b). Then constrained deformations are used to t the model (c). (d) shows the result obtained after having connected the di erent models through topological operations.

Fig. 14 .

 14 Fig. 14. Aneurism segmentation from an angiographic image.

  compares the result obtained by GCD (a) and B-spline registration (b). Since FFD has more DOF than B-spline registration, the GCD yields a more accurate t. This can be seen in sub gure (c) showing the distance to data. Moreover, GCD has a lower computational cost due to the prohibitive cost of B-spline computation (sub gure (d)).

Fig. 15 .

 15 Fig. 15. Comparison between the GCD and B-spline registration.

Fig. 16 .

 16 Fig. 16. (a) foot model (b) range data and deformed model (c) distance to data.

  Figure 17(b) shows the convergence of the algorithm for six di erent models during the rigid registration step.

Figure 17 (

 17 Figure 17(c) shows the convergence of the GCD step (with = 30%).

Fig. 17 .

 17 Fig. 17. Liver convergence measures.

  The registration computation cost is the sum of the cost for nding the closest points and the LSE algorithm: C registration = C closest (N) + C LSE (D) where N is the number of model vertices and D the number of transformation DOF. C closest (N) is only dependent on the numberof vertices. C LSE (D) depends on the number of DOF and the LSE algorithm used. Thus C closest (N) is constant for di erent transformations, and C registration is almost constant as long as C LSE (D) remains negligible compared to C closest (N) this is the case for rigid, similarities

Acknowledgment

We thank IRCAD and Dosigray for abdominal CT-scans, Brigham and Women's Hospital for brain MR images, and General Electric, Medical Systems Europe, for the 3D images of aneurisms (Courtesy of Dr. Yves Trousset).