H Delingette

M Hebert

K Ikeuchi

Image Segmentation and Shape Representation Using Deformable Surfaces 1

We present a technique for constructing shape representation from images using free-form deformable surfaces. An object is modelled as a closed surface that is deformed subject to attractive fields generated by input data points and features. Segmentation is achieved by initializing the surface at some location in the scene and by letting it deforms itself until it fits the input data. Surface deformation is controlled by applying standard theory of mechanical systems. The algorithm is general in that it makes few assumptions on the type of features, the nature of the data and the type of objects. The algorithm does not assume that a clean segmentation of the input data is available. It correctly recovers shape even in the presence of spurious features and data points. We present results in a wide range of applications: reconstruction of smooth isolated objects such as human faces, reconstruction of structured objects such as polyhedra, and segmentation of complex scenes with mutually occluding objects. We have tested the algorithm using data from different sensors including grey-coding and laser range finders and video cameras, using one or several images. We briefly describe the application of this representation to the problem of computing stable grasp position for the manipulation of unstructured objects.

Introduction

The recovery of object shape from 3D data is one of the key issue in vision. One could define this task as the segmentation of a large set of data points into shapes corresponding to objects in the scene. The shape representation should be general enough to handle a wide variety of scenes yet simple enough to be usable for other tasks such as recognition and manipulation. In other words, the shape representation should have enough parameters to describe the specificity of the shape but must have as few parameters as possible to be usable and to be robustly extracted from visual data. This conflict is related to the scale space problem, where one would like to find a description fine enough to capture the key details of the shape, but coarse enough to get rid of spurious details that would only penalize applications that use the description, and that would render the shape extraction process less robust.

In this paper, we propose an approach that attempts to resolve this conflict by substituting to fine / coarse opposition, the feature / data opposition. Several psychophysical experiments have proved that the human eye is able to capture the main shape of an object by seeing only a few characteristic elements or features. These features can be either geometric (distance discontinuities, surface orientation discontinuities, corners, minimum of curvature, etc.) or higher level such as reflectance properties. While these features capture most of the shape information, it is difficult without apriori knowledge to construct a full reconstruction of the object.

Several solutions have been proposed. Besl and Jain 6 built curvature-based object representations by classifying surfaces according to the sign of its principal curvatures. Pentland [START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF] presented a physically-based algorithm, to recover in a unique manner a model from a set of features and a set of vibration modes. In this approach, the key is to find the correct features, thus restraining the algorithm to either smooth or structured shapes. Another approach is to use both features and range data in separate stages. In a first stage, features are grouped into hierarchical sets, according to geometric properties (symmetry, connexity, etc.) and in the second stage models are fit to the segmented parts This idea of hierarchical representation was initiated by Marr and Nishihara [START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF] and pushed further by the seminal ACRONYM [START_REF] Broooks | Symbolic Reasoning Among 3D Models and 2D Images[END_REF] vision system by using generalized cylinders. Pentland [START_REF] Pentland | Recognition by Parts, Proc. ICCV[END_REF][START_REF] Pentland | Perceptual Organization and the Representation of Natural Form[END_REF] 's "Representation by parts" using deformed superquadrics, proved to have some successful results but encounters some limitations. While the feature grouping requires some accurate feature extraction and high level reasoning, the fitting of superquadrics [START_REF] Yokoya | Recovery of Superquadrics Primitives from Range images by Simulated Annealing[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF][START_REF] Pentland | Recognition by Parts, Proc. ICCV[END_REF] to the range data has some unstable behavior, due to its non-linear nature, and is suitable for only smooth and simple shapes.

Those techniques attempt to represent all shapes by using a set of elementary shapes (superquadrics, generalized cylinders, parametric patches, etc.) that can be described by a few parameters. This is clearly beneficial from the point of view of object recognition which amounts to manipulating analytical equations of the elementary shapes. In practice, however, it restricts considerably the class of objects and scenes to which the techniques can be applied. The shape descriptions will always be rough approximations unless the object being modeled can be exactly described by one of the elementary shapes. More general representations could be obtained by adding degrees of freedom to the elementary shapes (e.g., adding tapering and bending to superquadrics). However, the non-linear fitting algorithms involved in the recovery of such shapes become rapidly computationally expensive and numerically unstable. Furthermore, most of those techniques assume that the observed scene is first segmented in regions corresponding to individual objects. The shape extraction algorithms are then applied to each object. However, accurately segmenting a scene is a hard problem in itself. Therefore any shape representation algorithm that requires a perfect segmentation is doomed to fail in realistic situations.

To address those problems, Terzopoulos and Witkin 22, 16 , 22 proposed the concept of deformable contours and deformable surfaces. In this approach, contours or surfaces are deformed under the influence of forces generated by image elements such as edgels. Deformable contours are used for extracting 2-D shapes, deformable surfaces are used for extracting 3-D shapes from range and intensity images. In the latter case, for example, they used a tube made of elastic mate-rial that was deformed by a potential field computed from an intensity image. Using symmetry and smoothness as internal constraints, this tube is able to approximate the main shape of the object. But to perform correctly, the tube had to be initialized close to the line of symmetry because the deformations occur only locally. This work demonstrated that deformable shapes are powerful tools for image segmentation and shape representation.

Our approach is also to represent shapes by deformable surfaces. We consider surfaces that are topologically equivalent to the sphere. Given an observed scene, a surface is initialized in the vicinity of detected features. The surface is then deformed subject to forces generated by features and data points. The forces generated by data points control local shape, while forces generated by features control global shape. A smoothness energy is added to the deformation equations to take into account the fact that data and features may be sparse and noisy. Clearly, such freeform surfaces can represent a large class of objects since few constraints are put on the resulting shape. However, our aim is also to build a shape extraction algorithm that does not depend on the nature of features and data. Furthermore, the algorithm must be robust enough so that it does not require a precise segmentation of the scene as an input. The deformable surfaces algorithm presented in this paper satisfies those goals. Specifically, it has the following characteristics:

• Physically based algorithm: The representation is obtained by deforming a surface. The dynamics of the deformation is modelled by the Lagrangian equations of mechanical systems. The surface is constrained by external actions from features and range data and by internal actions. As demonstrated in previous works, using physically based procedures provides better control on the stability and convergence of the algorithm.

• Optimal shape description: Different pieces of information can be given different weights in the algorithm. For example, a strong feature like a corner would have more influence on surface deformation than data points measured in a featureless area. This leads to an optimal shape description in that it realizes the best compromise between different shape clues such as different types of features and data points with varying levels of noise.

• Stability: because the algorithm uses both features and data, it is less sensitive to spurious features, noisy data or missing data. Moreover, this stability enables to perform segmentation by discarding the features and data that is incompatible with the current shape of the surface. Therefore, our algorithm does not require a clean segmentation of the scene, an approximate partition of data and features into region of interest is sufficient. This is in sharp contrast with other techniques that require that the observed scene is already segmented into regions that correspond exactly to individual objects. Furthermore, the stability of the algorithm allows us to deal with non-uniform data distribution and heterogeneous features.

• Generality: The algorithm makes few assumptions on data and observed objects. The only requirement is that some features can be extracted from input data and that the minimum and maximum sizes of the object expected in a typical scene are known. We have applied successfully the algorithm to the reconstruction of smooth isolated objects such as human faces, to the reconstruction of structured objects such as polyhedra, and to the segmentation of complex scenes with mutually occluding objects. We have tested the algorithm using data from different sensors including grey-coding and laser range finders and video cameras, using one or several images. The generality of the algorithm allows us to deal with this array of applications without modifying its core. All that needs to be changed is the front-end program that converts feature and data from the sensor to an internal data structure. Section 2 describes the theory of deformable surfaces, Section 3 describes the implementation for various sensors and applications.

Deformable surfaces

The basic surface description is a free-form deformable surface. The surface is subject to a field of forces resulting in a total surface energy. The energy changes over time as the surface is deformed. A deformable surface is subject to three types of forces, external, internal, and inertial forces.

External forces are generated by input data points and input features. External forces apply deformations that bring the surface as close as possible to the data. External forces must be carefully designed to obtain an optimal shape description that uses features to constrain global shape and data points to fine-tune local shape ( Figure 2 ). External forces can be represented as links between the surface and data points and features ( Figure 2 ). Earlier work on surface reconstruction using splines subject to constraints from data points required a complete mapping between data and surface points to be defined beforehand (see Figure 1 ). This constraint restricts the geometry of the problem and requires the surface to be initially very close to the solution. To remove those constraints, we need to design the forces so that there is no restriction on their geometry.

Internal forces are generated by the surface itself as it is deformed. The inclusion of internal forces ensures that the surfaces will not tear apart, fold onto itself, or exhibit high curvature points or sharp discontinuities in curvature. The other role of internal forces is to provide constraints in regions in which little or no data is available. This is similar to the regularization approach to surface reconstruction from sparse data [START_REF] Sato | Range-Imaging System Utilizing Nematic Liquid Crystal Mask[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF] . The standard way of defining internal forces is to define the corresponding energy as the integral over the surface of the magnitude of the first and second derivatives [START_REF] Pentland | From Features to Solids[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF] which characterizes the surface smoothness. The relative importance of external and internal forces is a trade-off between accuracy and smoothness. High internal forces generate a very smooth surface that may be far from the input data. Low internal forces allow the surface to fit the data closely but they also allow the surface to fit any noise in the data. Furthermore, they do a poor job at interpolating regions with sparse data. Inertial forces are generated by the motion of the surface as it evolves over time assuming that the surface has a non-zero mass. Inertial forces are necessary to model the deformable surface as a dynamic mechanical system.

In this section we first describe the components of the three types of energy and the equations of motion (Section 2.1), then we describe in detail how each energy is computed (Sections 2.2, 2.3, 2.4) in the case of a continuous deformable surface. In Section 2.5 we describe the computation of the energies and the implementation of the equations of motion in the case of a discrete deformable surface.

General equations of motion.

The internal and external energies involved in the deformation of the surface are:

• Smoothness energy E smoothness : The smoothness energy is a measure of the average curvature of the surface. It is necessary to guarantee that the surface is reasonably smooth, especially in the case of sparse data. The smoothness energy is internal in that it depends only on the shape of the surface in the vicinity of each point.

• Feature energy E feature : The feature energy quantifies the effect of the features on the surface. Each surface point is attracted by each feature. The magnitude of the attractive field is function of the distance between surface point and feature.

• Data energy E data : The data energy quantifies the effect of data points on the surface. Theoretically, each surface point is subject to an attractive force by each data point. However, in practice the closest data point is the only one taken into consideration.

To calculate the equilibrium position of the surface using mechanical systems theory, we need to introduce two additional inertial energy terms:

• Kinetic energy T: A mass μ is associated with each data point thus generating a kinetic energy term. Unlike other dynamic splines [START_REF] Pentland | From Features to Solids[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF] , our scheme is explicitly linking the deformation of the surface with the minimization of the potential energy.

• Raleigh Dissipation energy D: The dissipation term is added to simulate the exchange of energy between the dynamic surface and a virtual medium in which it evolves. This damping term is added to avoid cases in which the surface oscillates forever around an equilibrium position.

Following the equations of mechanics and the principle of least action [START_REF] Yokoya | Recovery of Superquadrics Primitives from Range images by Simulated Annealing[END_REF] , the surface reaches a stable equilibrium when the Lagrangian of the system of forces reaches a minimum. The Lagrangian of the system is: L = T -D -E smoothness -E feature -E data . A deformable surface is parametrized by two spatial coordinates (u,v) and by the time t. A surface point P has coordinates: With these notations, using the calculus of variations on L(u,v,t), one gets: where μ is the mass density of the surface and k is the damping factor. This surface is being deformed between t=0 and t=T 0 . Thus, unlike other scheme using deformable splines [START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF][START_REF] Pentland | From Features to Solids[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF][START_REF] Pentland | From Features to Solids[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF][START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF] in which the surface is deformed until an error is less than a threshold, we control the amount of deformation of the surface by imposing a time limit T 0 . This allows us to deal with noisy and perturbing data.

MORE ON THE MYSTERIOUS T0

Smoothness Energy

We use the standard Tikonov's stabilizers of the first and second order as a smoothness measure. The first order measures the distance discontinuities while the second order measures the surface orientation discontinuities. Denoting partial derivatives by subscripts (e.g., , the energy is defined by:

r u v t , , (
)

X u v t , , ( ) Y u v t , , ( ) Z u v t , , ( ) = μ t 2 ∂ ∂ r u v t , , ( ) × k t ∂ ∂ r u v t , , ( ) × - F smoothness F feature F data + + ( ) + = r u u ∂ ∂ r u v t , , ( ) =
The corresponding force is therefore:

The coefficients α 1 and α 2 controls the relative smoothness of the surface.

Feature Energy

Because every feature contributes to the global deformation of the surface, our approach is to link every feature to every point of the surface. The links can be seen as springs attached to each feature. Feature energy is given at a surface point by the sum of the energies generated by all the springs attached to the features (see Figure 3 ): With this definition of feature forces, the algorithm satisfies two requirements:

• No initial mapping is needed, so that the initial surface can be "far" from the object.

• Due to the springs actions, the points are going to concentrate toward the features, thus yielding an optimal global description.

To get both global and local deformations is hard because the two types of deformation, feature and data, should be balanced. If only feature deformation were applied one would get a surface that connects the features with mostly planar surfaces in between because of the smoothness constraint. On the contrary, in the case of data deformation, one would get a surface that would faithfully follows the shape of the object only where the surface is close enough initially(see Figure 4 ). Our solution is to change the relative influence of both types of deformations over time: initially, the surface is mostly submitted to feature forces, putting the surface close to object; then the feature forces decrease to be replaced by data forces. This amounts to decreasing over time the distance within which the points of the surface are attracted by features. In the current implementation, we model features as 3-D line segments. We denote the midpoint of the segment by , and its length by . With these notations, the energy field generated by feature i, is:

where:

• U(x) is a function that is quadratic if x<1 and constant if x>1. U becomes constant after x = 1 to avoid large attractive forces generated by features that are very far away. ----------------------------

E smoothness α 1 r u 2 r v 2 + ( ) α 2 r uu 2 2 r uv 2 r vv 2 + + ( ) × + × [ ] u v d × d 0 1 ∫ 0 1 ∫ ⎝ ⎠ ⎜ ⎟ ⎛ ⎞ t d 0 T 0 ∫ = F smoothness α 1 Δ uv X Δ uv Y Δ uv Z α 2 Δ uv Δ uv X Δ uv Δ uv Y Δ uv Δ uv Z × + × = E feature E feature i i 0 = n ∑ = F i l i E feature i E feature i U l i PF i D 0 i t 0 t - ( ) × -
- ⎝ ⎠ ⎜ ⎟ ⎛ ⎞ u v d × d v ∫ u ∫ t d 0 T 0 ∫ =
This improves the robustness of the algorithm since spurious features that are not part of the object will not contribute significantly to the forces.

• is an approximation of the integral of the Euclidean distance between and the feature points over the segment, normalized by the length of the segment. Other feature representations can be used in the same framework by replacing this term by the appropriate value. For example, for a point feature it would be the distance between surface point and feature point.

•s is the distance between the feature and the center of the surface at its initial position. If the surface is initialized as a sphere, it is the distance to the center of the sphere.

Figure 5 shows how the potential field varies over time in the case of four coplanar features. The center of the surface coincides with the center of gravity of the features in this example. This example demonstrates that this dynamic potential avoids the "local minimum" problem, because there is no fixed valleys in which points can be trapped.

From the definition of feature energy, the force generated by the feature i at surface point P at time t is:

This force is roughly equivalent to the force of a spring if is less that and is null otherwise. Forces from each feature may be weighted to reflect the relative importance of different types of features.

Data Energy

Since data information is used only for local deformation, the corresponding force should decrease sharply with distance. Therefore, a gravity-type field instead of a spring model is appropriate for this type of deformation. Theoretically, a surface point is subject to forces from all the data points. However, for computational reasons we take into account only the closest data point. For every point P of the surface, the closest data point is denoted by P data . If W is the gravity potential, the resulting data energy is given by: Where K is a normalizing constant that has the dimension of a distance. Intuitively, surface points are influence only by data points that are at a distance less than K. Choosing K small enough compared to the maximum expected object size guarantees the locality of forces generated by data points. A gravity potential should be in 1/distance everywhere. To avoid infinite potential when the distance vanishes, the potential becomes a quadratic function near 0 (see Figure 7 ). The corresponding force is: ------------------- -------------------------------PP

l i PF i r D 0 i F feature i l i D 0 T 0 t - ( ) PF i × × ---------------------------------------------------- x d d × U l i PF i D i 0 T 0 t - ( ) × --------------------------------- ⎝ ⎠ ⎜ ⎟ ⎛ ⎞ ( ) PF i × = l i PF i D 0 i T 0 t - ( ) × E data W PP data K -
- ⎝ ⎠ ⎛ ⎞ u v d × d v ∫ u ∫ t d 0 T 0 ∫ = F data W' PP data K --------------------- ⎝ ⎠ ⎛ ⎞ PP data K --
data × =

Implementation

We have assumed so far that our model is a continuous surface topologically equivalent to a sphere parameterized in (u,v,t). In practice, however, we can manipulate only discrete surfaces. This raises the problem of the parameterization of such surfaces, and, in particular, the impossibility to map a sphere into a square in a uniform way. To avoid creating poles, we adopt the tessellated icosahedron as a structure. Each face of the icosahedron is subdivided to yield arbitrary resolution of the parameter space. The number of faces of the tessellation is , where N is the frequency of the subdivision. Typically, we use N=5 yielding a decomposition of the parameter space into 500 faces. Figure 8 shows the tessellation of the sphere for N=5. We use the center of each triangle as a node, every node having therefore three neighbors. The time axis is also discretized. The deformable surface is updated at unit time increments. If we write as the position of the node i at the time t, then the discrete version of the motion equation is:

The surface is initialized as a sphere at t=0. The surface is deformed by applying repeatedly the equation of motion at each node. Forces and are computed at each node independently in a parallel manner.

is computed by approximating the first and second derivatives of the surface by finite differences. The computation of can also be carried out at each node independently except that it involves neighboring nodes. The most expensive part of the algorithm is the computation of the closest data point P data used in the computation of . It is theoretically in O(lm) where l is the number of data points and m the number on nodes. In practice, some precomputations and assumption allow to improve greatly this computational time. The preprocessing depends on the nature of the data and is described in the next Section for each experiment. This preprocessing does not affect the deformable surfaces algorithm itself since it is part of the front end data preparation. The algorithm is otherwise linear in the number of features and the number of iterations. Due to its highly parallelizable nature, substantial speed-ups can be achieved.

Several parameters must be set to apply the motion equation. It is important to note that most parameters are independent of the nature of the data. The parameter settings in the current implementation are:

•Radius and center of initial sphere: The determination of center and radius of the initial sphere depend on the nature of the data. The initialization procedure is described in the next section for each type of data.

•Number of iterations, T 0 :

The number of iteration is T 0 (between 0 and T 0 -1) and is related to how far the initial shape is to the object; the largest T 0 is and the smallest the deformation is between two iterations. We currently use .

•Smoothness coefficient and : The smoothness coefficients should be between 0 and 1 with . Actual "good" values have to be determined empirically. We use and in our experiments.

•Damping factor k: The damping factor should be close to (but lower than) 1 to ensure smooth deformation of the surface over time and to avoid oscillations. We use .

•Normalizing factor K: K depends on the environment. It is computed as one-fifth of 20N 2

X t i X t 1 + i X t i 1 k - ( ) X t i X t 1 - i - ( ) × F smoothness F data F Feature i i 0 = n ∑ + + + + = F data F feature i F smoothness F smoothness F data T 0 200 = α 1 α 2 α 1 α 2 > α 1 0,75 = α 2 0,4 = k 0,9 =
the radius of the largest object expected in a scene for a given application. The algorithm is robust with respect to K so that a rough estimate of object size is sufficient.

Experimental Results

In this section we present experimental results obtained by applying the deformable surface algorithm to real image data. We present two sets of results. The first set is obtained using range data from an light-stripe range finder. Those experiments show that the algorithm can extract 3-D object shape from range images without requiring good feature extraction or perfect scene segmentation. The second set of experiments involves applying the algorithm to intensity images. The goal of those experiments is to validate the claim that the algorithm is independent of the nature of the data. We show that the algorithm correctly reconstructs shapes from intensity images.

Range Data

For our experiments with range data, we use a commercial light-stripe range-finder that consists of a camera and a projector that projects patterns commanded by a LCD board [START_REF] Zucker | The Organization of Curve Detection: Coarse Tangent Fields and Fine Spline Coverings[END_REF] . The sensor processes the images of the patterns using standard light-striping geometry. It delivers a set of four images of 240 rows by 256 columns, one intensity image and three images of the three coordinates of every pixel with respect to a reference frame. Several sensors can be used at once to yield multiple views of a scene. A calibration procedure is used to express all data points coordinates with respect to a single world-centered frame. We conducted the experiments with either single views or multiple views. Using multiple views demonstrates that the algorithm is completely independent of an image-centered reference frame. In particular, it does not use the grid structure of the image or the uniform sampling of pixels in the image. This shows that it can be used with non-imaging sensors that measure non-uniform sparse data. This is a major difference with other reconstruction algorithms that use the image grid as a discretization of the surface.

We analyze the experiments in a simple to complex fashion. We first discuss two experiments in which images of isolated objects, a cube and a human face, are used. Those two experiments are intended to show the basic operation of the algorithm. The choice of the two objects is not arbitrary. The example of the cube shows that the algorithm can recover the shape of objects with corners. At the other extreme, the example of the human face shows that the algorithm can recover the shape of a complex smooth object, including small local shape features. We then describe experiments in which scenes with multiple occluding objects are segmented and the object shapes are recovered. We use two examples. One is the segmentation of natural scenes that contain rocks and pebbles. The second one is the segmentation of a scene with three cubes stacked together. We use simple known shapes in the second example to be able to evaluate the quality of segmentation and shape recovery. In those two examples, no segmentation is applied prior to applying the deformable surface algorithm. Only a rough guess of where the objects are in the scene in needed. This demonstrates that the algorithm performs image performs scene segmentation as well as surface reconstruction.

Isolated Cube

Figure 9.a,b shows two range images of a cube that have been taken from opposite directions. Two types of features are computed, distance discontinuities and surface orientation discontinuities by using standard edge operators. The edge magnitude image is thresholded and a thinning algorithm is applied to the resulting regions. After polygonal approximation, the final feature representation is a set of 3-D segments. Figure 9.d shows the segments. Because of the geometric nature of the object, we are able here to extract correctly most of the edges of the cube.

The deformable surface is initialized as a sphere (see Figure 9.e); the sphere is set so that one diameter passes through a given point and is parallel to the line of sight of the camera (Figure 11). This procedure has the advantage of requiring only the coordinates of one point (the line of sight is given by a calibration file), thus allowing a fully automatic procedure if a high-level module can determine such point.

The result of the deformation is shown in Figure 9.f and Figure 10. Figure 10 shows the triangulated surface (every node has three neighbors). Figure 9.f gives a synthetic aspect of the reconstructed object using ray-tracing. It can be seen that the algorithm correctly recovered sharp edges and corners, even though they are slightly rounded due to the smoothness constraint. This is due to the optimal description capability of our algorithm: it differentiates information from features information from range data.Figure 10 shows clearly that those edges correspond to a concentration of nodes.

Human Face

A human face is a good example of a complex object with parts of very distinct nature: the forehead and jaws areas are of little interest for face recognition, whereas eyes, nose, mouth and chin are the main characteristics of a face. Building compact and accurate surface representation is crucial in the area of face recognition [START_REF] Cerrada | A 3D Object Reconstruction System Integrating Range Image Processing and Rapid Manufacturing[END_REF] . We would like a surface reconstruction algorithm to smooth the input data and generate a compact representation of the face while keeping an accurate description of the areas of interest. We applied our algorithm to a range image of a human face to verify that it does perform as desired. Figure 12.a shows the intensity image and Figure 12.c shows the range data of a face, the data corresponding to the hair of the person is removed because the sensor generates extremely noisy measurements there. Only a partial view of the face is available because of self-occlusions. For example, the left part of the nose is not visible. In the first experiment, we consider as features only the boundary of the face. Figure 12.c shows the segments that surround the face. Using this minimum information, the deformation of the surface initialized as a sphere (Figure 12.e) gives the shape shown in Figure 13.a. While the overall shape of the face was found, important facial features such as nose, mouth and eyebrow were smoothed so that only the general shape is recovered.

To avoid this smoothing effect it is necessary to tag the eye, nose and mouth as being important features. The thresholding of the magnitude of an edge detector on the intensity image provides a way to extract those features (Figure 12.b). A simple thinning process generates a new set of segment features that are used for a second experiment (Figure 12.d). The result of the final recovery is shown in Figure 13.b: this time the nose, eyebrow and mouth are clearly visible and the left part of the nose was interpolated. Several conclusions may be drawn from this experiment. First, it shows that fine details as well as general shape can be recovered due to the combination of feature and data forces. Second, it shows that the input features need not be very clean. Since a simple thresholding/thinning operator was used, no attempt was made here to extract clean edge features. Third, the algorithm provides a way to compress the input data drastically while retaining the important information. By using approximately 250 points (500 nodes for the whole surface) we are able to construct a fine and accurate description of the face.

Stacked Pebbles

The generality of our approach makes it particularly well-suited for unstructured environment in which there are few constraints on object shapes. The scene we use in this example consist in a set of stone pebbles that are lying on top of sand. In order to have more reliable data and features, we took three range images from three different viewpoints. The features are of three natures: distance discontinuities, surface orientation discontinuities and shadow boundaries. Shadows are parts of the scene that are visible from the camera but are not illuminated by the projector. Shadows are used here as clues for the presence of an object since a shadow is created by an object occluding the light coming from the projector. The feature segments are shown in Figure 14.b. The figure demonstrates that using classical feature-based segmentation techniques would be a very challenging task.

The segmentation proceeds by selecting a shadow region, starting with the largest in the image. Using the geometry of the sensors as computed from an off-line calibration procedure, a 3-D point that is assumed to be inside the object is selected based on the position of the shadow in the image. The center of the initial sphere is initialized at this point, its radius is initialized to a constant value that is the average size of the expected objects. Figure 14.c shows the initial sphere. Data points and features that are less than a fixed distance away from the starting point are taken into account in the computation of the forces. The distance is computed from the radius of the largest object expected in the scene. Figure 14.d,e show the shape as it evolves. Figure 14.f shows the final shape. This experiment validates several of our claims. First, shape is recovered even though there are few apriori constraints on object shape. Second, weak feature extraction is sufficient to recover shape. Most importantly, no scene segmentation is required. All that is needed is to hypothesize one point inside each object based on the knowledge of shadow geometry for this particular sensor. No feature-based segmentation is necessary. The algorithm is able to deal with numerous features and data points that are not part of the object. Finally, by using multiple views, this experiment shows that the algorithm works directly in 3-D and is independent of the image structure.

This application to the segmentation of natural scenes is part of an operational system that automatically picks up pebbles in cluttered natural environment [START_REF] Francois | A Three-Finger Gripper for Manipulation in Unstructured Environments[END_REF] . The discrete surface extracted from range images are fed to a program that computes the stable grasp positions for a three-finger gripper. A manipulator executes the grasping. Figure 15 shows the gripper picking up a pebble in a typical environment. This shows the relevance of the algorithm in applications in which accurate surface descriptions of arbitrary objects are needed. In particular, the grasp determination algorithm uses the description of the surface to compute local properties such as curvatures. This would not be possible using standard parametric description such as superquadrics.

Stacked Cubes

To investigate further the ability of the algorithm to perform segmentation and surface reconstruction, we designed a simple scene with three overlapping cubes. Using simple objects of known shapes such as cubes allows us to compare the resulting segmentation and reconstruction with the actual scene more easily. Figure 16.a,b shows two opposite views of a scene where one cube is in balance between two adjacent cubes. The shape recovery was applied to each of the three cubes. The square drawn in the Figure 16.a, b, corresponds to the region that is used for the recovery of the first cube in both images. All features and data points that fall within this region of influence are taken into account in the recovery process. Therefore, the objects are recovered without explicit segmentation of the images. As before, the surface is initialized as a sphere of fixed radius whose center is roughly at the center of a group of features. The center of the sphere is also the center of the region of influence. Figure 16.c to Figure 18.a show features, range data, initial surface position, and a shaded display of the recovered surface for each of the three objects. Figure 18.b shows the final description of the scene displayed as three shaded surfaces. This experiment shows that the scene was correctly segmented into objects with the correct shapes even though no explicit segmentation was done except for the selection of region of influence. Among the three reconstructed shapes, the second seems to be the best (Figure 17.d), mostly because he is not occluded by any other cube. On the contrary, the third cube (Figure 18.a) has one face that is poorly recovered. A look at the range data shows that there is no feature and no data corresponding to this face of the cube, most probably because of the occlusion from other cubes. In this case, the algorithm does the best job it can using only the smoothness energy. The reconstruction of the scene clearly demonstrates that all visible parts were correctly extracted. By computing for every node the minimum distance with range data, it is possible to identify which point was interpolated or not, and therefore to use only the reliable parts of the reconstructed surface later on.

Conclusion

Figure 1 .Figure 2 Figure 3 FeatureFigure 4 a,Figure 6 a,

 12346 Figure 1 . Nodes of a spline are attached to data points via springs

Figure 9 a,

 9 Figure 9 a, Intensity image with features from the first view; b, Intensity image with features for the second view; c,Range data from thetwo views combined; d, Features from the two views combined; e, Cut-section of the initial position of the surface; f, Final shape of the surface;

Figure 11 Figure 13 a,Figure 14 a,Figure 16 a,Figure 17 aFigure 18 a,

 111314161718 Figure 11 Initialization of the sphere

29 29 29 29 29 29 29

1. This research was supported in part by NASA under Grant NAGW 1175, and in part by DARPA through ARPA order No. 4976. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied of NASA, DAR-PA, or the US government.


Radu Stefanescu

Xavier Pennec

Nicholas Ayache
email: ayache@sophia.inria.fr

Xavier Stefanescu

Nicholas Pennec

A Grid Service for the Interactive Use of a Parallel Non-Rigid Registration Algorithm of Medical Images

Keywords: grid service, non-rigid registration, cluster, medical imaging

Objective: The goal of this work is to improve the usability of a non-rigid registration software for medical images.

Method:

In order to use the interactivity of a visualization workstation and the computing power of a cluster, we have build a registration grid service. On the user side, the system is composed of a graphical interface that interacts in a complex and fluid manner with the registration software running on a remote cluster.

Conclusion:

Although the transmission of images back and forth between the computer running the user interface and the cluster running the registration service adds to the total registration time, it provides a user-friendly way of using the registration software without heavy infrastructure investments in hospitals. The system exhibits good performances even if the user is connected to the grid service through a low throughput network such as a wireless network interface or ADSL.

Introduction

Non-rigid registration recovers the deformation between two images. As an example, Figures (1a) and (1b) present two-dimensional slices of three-dimensional magnetic resonance images (MRI) of the brain coming from different subjects. The algorithm recovers a deformation, presented in Fig. (1c), that morphs one of the images (Fig. 1a) into an image (Fig. 1d) that is as close as possible to the other (Fig. 1b). Since the two brains are not identical, this correspondence can not be perfect. Moreover, it is necessary that the recovered deformation preserve the anatomic coherence of the objects in the images: the deformation has to be smooth and invertible. This non-rigid registration procedure is necessary to compare and/or fuse images in many domains of medical imaging, ranging from building and usage of medical atlases [START_REF] Dlg | A Dynamic Brain Atlas[END_REF] to pre-operative planning [START_REF] Warfield | Real-Time Registration of Volumetric Brain MRI by Biomechanical Simulation of Deformation during Image Guided Neurosurgery[END_REF]. Generally, non-rigid registration is a computationally intensive, and hence time-consuming procedure.

In order for registration to be really used in clinical environments, physicians need software tools that are precise and robust. Furthermore, care must be taken regarding the usability of such tools. One primary issue is the execution time, which should be lower than a few minutes. However, since financial pressure makes it easier for healthcare organizations to invest in medical equipment than in high-performance computing, the software should not require hardware that is too expensive to be acquired by a hospital (more expensive than medical equipment), or that is too hard to administrate.

In a previous work [START_REF] Stefanescu | Grid Enabled Non-Rigid Registration with a Dense Transformation and A Priori Information[END_REF], we developed a non-rigid registration algorithm able to reliably recover large deformations and that admits an efficient parallel implementation on a cluster of personal computers. The algorithm solves a two-part partial differential equation: if we denote by U=(u 1 ,u 2 ,u 3 ) the deformation that links the source image J to the target image I , then the first term minimizes the square distance

2 )) ( ( ) ( p U p J p I + -
between the target image and the deformed source image. The second term is based on two assumptions: the deformation is supposed to be regular, and the degree of regularity spatially varies according to the deformability of the imaged objects. Hence, the regularity term has the form of a diffusion term, where the conductivity d represents the local ability to deform. The two terms are weighted by a confidence factor k(p) that estimates at each point the reliability of the matching process:

[ ] ∑ Ω ∈ ∇ ∇ - + + - ∇ = ∂ ∂ ) ( 2 ) ( ) 1 ( )) ( ( ) ( images p i i u d k p U p J p I k t u
The parallel implementation of this algorithm, using the Message Passing Interface (MPI) and a cluster of personal computers, lowers the registration process time for typical 256x256x124 MR images to less than five minutes on a cluster of 15 2GHz Pentium IV PC s linked by a 1GB/s Ethernet network.

Though accurate, robust and fast, our registration software suffers from several issues regarding its usability. First, it is command line based, which means that once the algorithm is started on the two input images, there is no mean to interact with the algorithm until it ends. Desirable interactions are: stopping a registration process either because it is about to diverge or because the level of convergence is sufficient for the user s needs; modifying parameters; visualizing the intermediate result. In the near future, we would also like to enable the user to add some medical expertise in order to guide the registration process. All this requires the creation of a graphical user interface (GUI) that allows the user to master the registration software.

Our system is composed of two interacting modules: the registration software and the graphics user interface. The registration module performs the computation. It is run on a parallel computer that usually needs special conditions to operate, such as air-conditioned rooms. Therefore, such a machine cannot lie next to the user and is more likely to find a suitable place as a shared resource in a data center than in an operating or pre-operative planning room. The GUI provides the interaction between the registration module and the user. It also deals with the visualization of the three dimensional images involved. For visualization performance reasons, the GUI must be executed on the user s computer. The computers running the two components of our system must be connected through a network. Due to the different conditions in which these two computers must operate, such a network may be a long distance one. Moreover, it may be difficult to install additional wires in the operating room for intra-operative registration. Thus, we want to enable our service to run through a wireless, hence low bandwidth, network interface.

The purpose of our work is to implement the registration as a grid service able to connect the clinically useful user interface to a computational center. This interaction has to be sufficiently refined to allow the user to dynamically manipulate the state of the registration software.

A previous work on a non-rigid registration service was done by Ino et al. [START_REF] Ino | A High Performance Computing Service over the Internet for Nonrigid Image Registration[END_REF], who used a high-end cluster connected to the visualization workstation through a high-bandwidth Wide Area Network (WAN). We believe that the infrastructure these authors used is not representative of the one available in most hospitals. Lower end to average clusters and low bandwidth WAN s are in our opinion more realistic hardware platforms for any medical grid service. Moreover, this paper does not address the interactive usability issue that we believe essential in a clinical environment.

In this paper, we present a system that enables the interactive use of parallel registration software running on a remotely located cluster. The presentation emphasizes the mechanisms used to connect our user interface to the grid service, and the security and usability concerns that are raised by such a system.

Objectives

To summarize, we want to build a registration grid service implemented on a cluster of PC s on one side, and a graphical interface able to control it on the user side. In order to connect them, we need a communication library that fulfills several conditions:

1)

The communication library has to be as flexible as possible. If the users makes two requests A and B to the GUI, and request A is transmitted to the grid service before request B, it is important that no constraint be imposed to the responses of the two request. The response to B may arrive before or after the response to A, and one or both responses may be never sent at all.

2)

Communications have to pass through firewalls.

3)

The access to the grid service should be secured. A public key authentication such as RSA fulfills our needs regarding access permissions to the grid service.

4)

Communications should be encrypted. Since confidential medical data is transmitted through a WAN that is not always under the hospital s full control, this data should be encrypted.

5)

The whole system should preserve the anonymity constraints imposed by current regulations [START_REF] Herveg | Directive 95/46 and the Use of Grid Technologies in the Healthcare Sector: Selected Legal Issues[END_REF]. Image files usually contain information about the subject s identity. This information should never make its way to the non-controlled WAN. The solution is to send only the data necessary for the registration: the image sizes and the image intensity data.

6)

In order to maximize the accessibility to the grid service, the communication library should spare bandwidth. Therefore, each message transmitted through the WAN should be compressed before sending.

Technological choice

During the recent years, a vast amount of work has been invested in distributed computing environments providing standardized and secure communications between remotely located computers. We review below some of these environments in view of our application.

Client/server

Several standardized communication methods between a client and a server have emerged. Recently, older standards, such as RPC and CORBA, have been replaced by protocols based on XML, which ensure the interoperability between implementations. The SOAP protocol, used by Web Services (WS) and the Open Grid Services Infrastructure (OGSI) [START_REF] Tuecke | Open Grid Services Infrastructure (OGSI) Version 1.0[END_REF][START_REF] Foster | The Physiology of the Grid: An Open Grid Services Architecture for Distributed Systems Integration[END_REF] has lately become the communication standard for client/server distributed systems.

Generally, the server interface is described in a dedicated interface description language (IDL in CORBA, WSDL in Web Services, GWSDL in OGSI), and the communication interface is generated automatically. The main advantage of such a system is that, when the client sends a service request to the server, the identification of the request and the decoding of its parameters is done in a transparent manner. We feel that several aspects of client/server systems makes them unsuited for our task: First, a client cannot simultaneously communicate with more than one server. This means that when the server is a parallel software, the client would have to communicate exclusively with the master node. However, when the master node receives a request, the identification of the request is done in a transparent manner, outside the user s control. Thus, further action would have to be taken in order to forward the request type to the other nodes. From this point of view, client/server would not bring additional functionality to the system.

The client/server system has another major drawback in our case. In this model, the server s sole purpose is to answer requests from the client. It cannot send requests to the client, and it remains idle between processing two successive requests. In our system, this constraint imposed on the server is too hard.

Message passing

This model, used by basic networking protocols (TCP/IP) and low level grid layers (Globus I/O), and generally adopted by computation libraries (MPI, PVM), provides flexible high-performance communications. In TCP/IP, the user has to explicitly describe the data transfers and the message encoding. Higher level libraries (MPI, PVM) facilitate the message encoding and also provide collective communications (e.g. broadcast). However, these libraries were designed for non-secured communication between the nodes of a parallel computer. Furthermore, they are generally unable to pass through firewalls.

We prefer the message passing communication model, since it does not impose any constraints on the possible actions of the service. Our goal is to create a message passing library that can provide secured communications and is compatible with firewalls.

Method

We designed a message passing library that is able to transmit messages back and forth between the GUI and each of the nodes running the grid service. Thanks to compression and encryption, the transmission of messages is fast and secure.

Communications

In our program, we distinguish between computation messages (exclusively exchanged between the nodes of the cluster) and control messages (exchanged between the parallel program and the GUI). Computation messages are directly communicated using MPI s pointto-point or collective communication subroutines, and they do not concern this article. Control messages are programmed as C++ classes containing a message tag that identifies the message, and a message handler describing the action to execute upon reception. Each message has two methods, pack and unpack, that describe the way the message can be packed into a buffer for sending, or unpacked from a received buffer.

Control messages are sent through communication channels, that encapsulate the communication protocols. They are C++ classes that provide:

• packing and unpacking basic data types (integer and real numbers and arrays) into buffers.

These functions are called by the message s pack and unpack procedures upon the sending and the receipt of the message. • sending and receiving the packed buffer.

Until now, we have implemented two channels:

• The long distance channel, currently based on TCP/IP, is used for communicating between the GUI and the master node. For performance reasons, this channel transparently compresses the messages before sending. It also takes care of the endianism difference that might exist between the machines running the GUI and the cluster. • The local channel, based on MPI, is used to forward the messages received by the master node from the GUI to the other nodes of the cluster.

The communication of a message using a channel is done as follows: the message tag and the message itself are packed into a buffer of the appropriate size and then sent. Upon receipt, the message tag is decoded and an object of the appropriate message class is instantiated. The message is unpacked and its associated message handler is invoked.

In the case where the GUI sends a message to the cluster, the message is first sent to the master node using the long distance channel, and then distributed by the master node to the entire cluster using the local channel (see Figure 2). The message handler is subsequently invoked in parallel by all the nodes of the cluster. The sending of a message by the service to the GUI is performed by the master node of the cluster.

The architecture we have presented here has the advantage of being modular. Adding message types or communication channel can be easily done by inheriting the appropriate C++ class.

Security issues

Both the GUI and the grid service are isolated from the rest of the world by firewalls. Regarding the long distance channel, three security issues are important:

• Communications must be able to pass through the two firewalls.

• The user must be authenticated before being allowed to use the service.

• Communications should be encrypted.

We use ssh tunneling in order to fulfill these three requirements. If a ssh tunnel is not created the registration service is not visible by the user. In order to create the ssh tunnel, the user has to authenticate using a user name and a password delivered by the service administrator. The tunnel also provides the encryption of the data flow between the GUI and the service.

Despite all these security mechanisms, a possible anonymity breach can come from the facial reconstruction that can be done on some MRI images. The solution would be to modify the images before transmission through the network in a way that does not affect registration but makes facial reconstruction impossible. We have not tackled this problem.

Results

The final system is the following: After user authentication, the GUI enables the loading and comparative visualization of the images to register. A user interface button enables the user to contact the grid service and start the registration. During the process, the user is informed in real time about the status of the algorithm (small size messages). An intermediate result image (several megabytes of data) is regularly sent to the GUI which displays it. At any time the user has the option to abort the registration. Upon termination of the algorithm, the registration result is sent back to the GUI.

We have tested our system by running the GUI on a 600MHz Pentium III laptop equipped with a wireless network interface. The grid service runs on 15 2GHz Pentium IV PC s linked through a 1GB/s Ethernet network and protected by a firewall running on a 2GHz Pentium IV PC. All systems run the Linux operating system. The resolution of the input images is 256x256x124, and their size is 47.2MB. After compression, they reach a size of only 3.3MB. The transmission of the input images takes 39 seconds, including the additional compression and decompression times. For comparison, if compression is not used, the transmission takes 2 minutes and 26 seconds. During the execution of the algorithm, each update of the intermediate result takes 18 seconds. The extra time took by the transmission of images through the network is of about 2 minutes, thus bringing the total registration time to about 7 minutes. However, thanks to the user interface and the grid service, the registration software becomes usable in the clinical environment. Thus, we believe that the two extra minutes in computation time are tolerable.

We also tested our system in a worst case scenario by using a home Asymmetric DSL connection between the user interface and the grid service. The maximum transfer rate is 16kB/s from the GUI to the service and 64kB/s from the service to the GUI. The image upload time is about 3min30, and the result update time is about 1min30. In this extreme case, the total registration time reaches ten minutes, including image transfer time. However, this proves that our non-rigid registration service can be used from almost everywhere with a very modest equipment.

Conclusion and future work

We have presented a system that enables a clinician to transparently use a computationally intensive but powerful non-rigid registration algorithm run by a parallel computer physically located at a large distance. The system combines the speed and precision provided by the parallel computer to the ease of use of the physician s workstation. A graphics interface enables the user to supervise the execution of the algorithm in real time. In the future, we will modify the system so that physicians will be able to use their clinical expertise to guide the registration. We will add an integrated authentication, the automatic creation of the ssh tunnel and the remote execution of the service by the user. A test in a clinical setup is also in preparation.

Our non-rigid registration service opens up numerous advanced medical image analysis applications to the clinical practice, such as the usage of brain atlases or image guided therapy. However, while medical imaging algorithmic issues have mostly been dealt with, grid management problems still have to be addressed. The most important problem is the standardization. Indeed, while the long distance TCP/IP channel fulfills its functionality requirements, its integration with standard security and authentication mechanisms and discovery services is not ensured. These aspects have been addressed in existing grid environments, such as the Globus Toolkit. In the future, we intend to replace the TCP/IP-based long distance channel with one based on Globus I/O. This would also enable us to replace SSL certificates and SSH tunneling with the similar but more standardized functionality provided by the Globus Security Infrastructure.

Figure 1 :

 1 Figure 1: The algorithm recovers the deformation (c) that morphs the 3D source image (a) into the 3D target image (b). The morphing result (d) should be as close a possible to the target (b).

Figure 2 :

 2 Figure 2: A message sent by the GUI is compressed, encrypted and sent through the ssh tunnel to the master node, which broadcasts it to all nodes of the parallel computer running the registration service.