Isaac Cohen
email: isaac.cohen@inria.fr

Isabelle Herlin
email: isabelle.herlin@inria.fr

Optical Flow and Phase Portrait Methods for Environmental Satellite Image Sequences

Keywords: Optical flow, Non quadratic regularization, Finite element method, Adaptive mesh, Phase portrait, Flow pattern classification, Ocean circulation

We present in this paper a motion computation and interpretation framework for oceanographic satellite images. This framework is based on the use of a non quadratic regularization technique in optical flow computation that preserves flow discontinuities. We also show that using an appropriate tessellation of the image according to an estimate of the motion field can improve optical flow accuracy and yields more reliable flows. This method defines a non uniform multiresolution scheme that refines mesh resolution only in the neighborhood of moving structures. The second part of the paper deals with the interpretation of the obtained displacement field. We use a phase portrait model with a new formulation of the approximation of an oriented flow field. This allows us to consider arbitrary polynomial phase portrait models for characterizing salient flow features. This new framework is used for processing oceanographic and atmospheric image sequences and presents an alternative to the very complex physical modelling techniques.

Introduction

Oceanographic and atmospheric images obtained from environmental satellite platforms present a new challenge for geosciences and computer vision. The wide ranges of remote sensors allow to characterize natural phenomena through different physical measurements. For example Sea Surface Temperature (SST), Altimetry and Ocean color can be used simultaneously for characterizing vortex structures in the ocean. A major advantage of environmental remote sensing is the regular sampling of the measurements and their availability. These regular temporal and spatial data samplings allow to characterize the short range evolution of atmospheric and oceanographic processes with image sequence processing.

In this paper we focus on dynamic satellite image analysis. The purpose of this paper is to derive a complete framework for processing large oceanographic and atmospheric image sequences in order to detect global displacements (oceanographic streams, clouds motion, ...), or to localize particular structures like vortices and fronts. These characterizations will help in initializing particular processes in a global monitoring system.

Processing such an image sequence raise some specific problems. Indeed, computing an apparent motion field to characterize short range evolution must take into account discontinuities of the motion field that occur near SST temperature fronts and clouds' boundaries. For this purpose we make use of a new regularization method for solving the optical flow constraint equation which involves a non quadratic regularization allowing flow discontinuities. Furthermore, an important issue is the reduction of the numerical complexity due to the large size of the images (typically 2048 x 1024). The minimization model is handled through a finite element method allowing the use of a non uniform domain tessellation. This tessellation is obtained with a mesh subdivision model allowing to obtain, locally, a finer mesh resolution near moving structures. These structures are detected and localized by studying the norm of the estimated motion field along the direction of the image gradient. Using such a non uniform domain tessellation yields an accurate optical flow near moving structures and a lower numerical complexity.

The second part of the paper deals with the interpretation of the obtained displacement field. A velocity field can be studied from two different view points. First, it is a motion field and therefore it can be studied according to the analysis of fluid motion. On the other hand, it is a vector field whose topology is described by its critical points and its salient features. This approach is more suitable in our context since characterizing salient flow features will help us in locating interesting structures (like vortices and fronts) which represent physical phenomena appearing in SST and atmospheric image sequences. We present in this paper a new approach for approximating an orientation field and characterizing the stationary points of the trajectories obtained from an arbitrary polynomial phase portrait. Furthermore the model is always linear, independently of the polynomial representation.

2

Non Quadratic Optical Flow Computation

The differential techniques used for the computation of optical flow are based on the image flow constraint equation:

dI(x, y, t) -I=u + I~v + I~ =0,

(1) dt where the subscripts x, y and t represent the partial derivatives. This equation, based on the assumption that the image grey level remains constant, relates the temporal and spatial changes of the image grey level I(x, y, t) at a point (x, y) in the image to the velocity (u, v) at that point [START_REF] Horn | Determining optical flow[END_REF]. Equation (1) is not sufficient for computing the image velocity (u, v) at each point since the two velocity components are constrained by only one equation; this is the aperture problem. Therefore, most of the techniques use a regularity constraint that restrains the space of admissible solutions of equation (1) ([1] and references therein). This regularity constraint is generally quadratic and therefore enforces the optical flow field to be continuous and smooth. But, true discontinuities can occur in the optical flow and they are generally located on the boundary between two surfaces representing two objects with different movements. This type of discontinuity occurs for example on temperature front in SST images, and cloud boundary in atmospheric images. Recovering this discontinuity is necessary for further analysis of oceanographic and atmospheric images. Indeed, locating and tracking of temperature fronts in oceanographic images represent an accurate estimation of the oceanic surface circulation. These fronts are defined as regions where the temperature variation is high. An accurate computation of optical flow components near these regions must take into account the flow field discontinuity along the temperature front. For this purpose, we defined a non quadratic regularization scheme preserving flow discontinuities while insuring a unique solution of equation (1). Several authors [START_REF] Black | Recursive non-linear estimation of discontinuous flow fields[END_REF][START_REF] Black | Robust dynamic motion estimation over time[END_REF] proposed non-quadratic schemes for the motion field regularity constraint. These schemes are based on the Graduated Non Convexity method, and allow a computation of optical flow components that preserves discontinuity but can not handle non uniform tessellation of the image which is an important issue when dealing with very large image sequences.

The proposed method makes use of the L 1 norm (defined by lull = /lul) for the regularization constraint. The advantage of this norm is that the variation of expressions like lull produces singular distributions as coefficients (e.g. Dirac functions). This property allows to preserve sharp signals as well as discontinuities in the space of L 1 functions. Such a property can be used to constrain the set of admissible solutions of Eq. (1). Considering the space of functions with bounded variation, i.e.: BV~ = {f =(fl,f2) such that ~ IVfll+lVf21dxdy< +oc}, the optical flow problem can be stated as the minimization of the functional :

~ 2 --~'~ (I~u 5) 2 dxdy,
where us and Uy (resp. v~ and v~) represent the partial derivatives of u (resp. v) with respect to x and y. This minimization problem (2) can also be viewed as a constrained minimization problem where we search for a vector flow field in the space of functions with bounded variations BV1, with the constraint of satisfying the optical flow stationarity equation (1).

The solution of the minimization problem ( 2) is obtained through the Euler -Lagrange equations:

{~u + (u*~ + vI=X~ + I=I,) = 0 Dv + (~z,I~ + vi~ + 5I,) = 0 ( 3
)
where D is the nonlinear operator defined by: ox Equation ( 3) is nonlinear and therefore must be processed in a particular way. An efficient method for solving this kind of nonlinear partial differential equations is to consider the associated evolution equations, or equivalently, the gradient descent method [START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF][START_REF] Rosen | The gradient projection method for nonlinear programming. Part I: Linear constraints[END_REF]. This time-dependent approach means that we solve the evolution equation:

{O~t +Du+(uI2+vI~I~+l=lt)=O + Vv + (uIdu + vI~ + &,) = 0 (5)
and a stationary solution characterizes a solution of Eq. ( 3). This evolution equation is solved through a finite element method allowing the use of arbitrary tessellations of the image domain by taking into account an index characterizing image motion in order to reduce the numerical complexity of the algorithm and increase its accuracy near moving structures (see Section 3).

3

Mesh Subdivision Computing an optical flow field over an image sequence using a classical approach leads to the solution of a large set of linear equations in the case of the quadratic regularizer or to an iterative solution for the non quadratic case. In both approaches space discretization (i.e. image tessellation) is an important issue since it defines the accuracy of the solution and the numerical complexity of the algorithm. In this section we propose a selective multi-resolution approach. This method defines a new approach for coarse to fine grid generation. It allows to increment locally the resolution of the grid according to the studied problem. The advantage of such a method is its lowest numerical complexity and its higher accuracy. Each added node refines the grid in a region of interest and increases the numerical accuracy of the solved problem [START_REF] Ciarlet | The finite element methods for elliptic problems[END_REF].

The method is based on a recursive subdivision of a given tessellation of the domain. This tessellation must fulfill the conform triangulation requirement of the FEM scheme [START_REF] Ciarlet | The finite element methods for elliptic problems[END_REF], i.e.: any face of any n-simplex 7"1 in the triangulation is either a subset of the boundary of the domain, or a face of another n-simplex 7"2 in the triangulation. This requirement restrains the type of n-simplex that can be used for an automatic non-uniform cell subdivision. We consider only triangular cells which are well adapted for domain triangulation and allow to derive a simple recursive subdivision scheme. Given a triangulation T = [.Ji T/ of the image domain, we have to refine it only near moving structures. These structures are characterized through w• representing the norm of the estimated motion field along the direction of the image gradient. Such a characterization is based on the image spatio-temporal gradients and does not require the knowledge of the optical flow field.

Let I represents the image brightness and VI its gradient vector field, then the optical flow equation (1) can be rewritten as [START_REF] Verri | Motion field and optical flow: Qualitative properties[END_REF] OI/Ot VI w--Y= IIVlII liVlll [START_REF] Ford | Nonlinear phase portrait models for oriented textures[END_REF] representing the component of the motion field along the direction of the gradient in term of the partial derivatives of I. Although w• does not always characterize image motion due to the aperture problem, it allows to locate moving points. Indeed, w• is high near moving points and becomes null near stationary points. The definition of w• gives the theoretical proof of the motion measure D defined by Irani et al [START_REF] Irani | Detecting and tracking multiple moving objects using temporal integration[END_REF] and used by several authors [START_REF] Muller | Adaptive-complexity registration of images[END_REF][START_REF] Szeliski | Motion estimation with quadtree splines[END_REF]: where W is a small neighborhood and C a constant to avoid numerical instabilities. This motion measure, defined as residual motion, is a particular form of w• where the numerator and the denominator are summed over a small neighborhood.

D(x,
The subdivision scheme is based on a splitting strategy. We start with a coarse tessellation of the image and split each cell T of the triangulation according to the norm w• summed all over the cell: w~. A cell is subdivided while w~ is greater than a given threshold and while its area is greater than another threshold. Furthermore, each subdivided cell must satisfy the conform triangulation requirement.

Experimental Results

We have experimented our method on the "Hamburg Taxi Sequence". In this street scene, there were four moving objects at different speeds. We first construct an adaptive mesh from the spatio-temporal gradient according to the algorithm described in section 3. We set the motion threshold, i.e. w~, to one pixel/frame and the cell area threshold to 15 pixels. The obtained mesh and the optical flow field are displayed in figure 1. We can observe that the mesh resolution is finer near moving objects allowing an important reduction of the algorithmic complexity, since we deal with 499 nodes and 956 triangles yielding an optical flow numerical accuracy lower than 0.5 pixel near moving structures. To achieve such an accuracy with a classical rectangular mesh, one has to consider 2560 nodes. The complete processing (i.e. the image gradients and the solution of Eq. ( 5)) takes 20 seconds on a Alpha 3000/500 workstation. Our main objective is to derive the surface ocean circulation from a sequence of Sea Surface Temperature measurements. These infra-red measurements are the most reliable for surface motion estimation and are daily available through the NOAA satellites. These measurements are corrupted by the presence of clouds, consequently, the data usually considered are composite images obtained by considering at each image pixel the maximal value during a period of time. Figure 2 shows a SST image which concerns the confluence region near Argentina coasts where a combined stream is formed by the Falkland northward flowing current (cold water), and the southward flowing Brazilian current (hot water). The discontinuities are located in the regions where these two streams are combined (top of figure 2). The obtained displacement field and mesh used are also displayed. ; ::::: : ::.i: .. ..; . . :'. \ "~ \ "~, X ><["q",, ~ \ ~\ \ \ "-q",,q\ ",-,,

9 t . p / / ~- # 1 1 . ~ ~ ~,
: [ . . . . : : : . : . . 9 Fig. 2. An illustration o~ the adaptive mesh approach for computing optical flow components from an image sequence. This figure displays a frame of the SST image sequence of confluence region near Argentina coasts. In this region discontinuities are located in the regions where the streams are combined.

E s t i m a t i o n o f a P h a s e P o r t r a i t

The computation of optical flow yields a quantitative measure of the flow field on each image point 9 By processing sea surface temperature and atmospheric image sequences we are also interested by the nature of motion since some phenomena like vortices are characterized through some specific patterns of the motion field 9

The extraction of higher level descriptors from a flow field is naturally crucial when studying fluid motion and vector analysis.

A classical approach is to use qualitative differential equations to characterize the orientation flow field by considering it as the velocity field of a particle in a dynamic system. Consider a particle governed by the model:

g(x,y)={~ =P(x'y) (8) = Q ( x , y)
This criterion is easier to handle since it is quadratic. In [START_REF] Cohen | A motion computation and interpretation framework for oceanographic satellite images[END_REF] we have proved that recovering the coefficients of the two polynomials P and Q E Q,~(~2) by minimizing $2 amounts to an eigenvalue problem.

This new formulation of the phase portrait similarity measure allows to derive a linear algorithm characterizing arbitrary polynomial portraits.

6

Flow Pattern Classification Analyzing a flow pattern consists in deriving a symbolic description from the model ( 8) fitted to the given local orientation field. In the case of optical flow, it would give a displacement information while in the case of fluid dynamic it will characterize the fluid velocity and the coherent structures of the flow. For example, the dynamic of the clouds can be tracked with infrared measurements, one can use such an approach in order to obtain a qualitative description of the atmospheric circulation without using a complex physical model of the underlying phenomena (see Figure 3). The symbolic description of the trajectory of the phase portrait model was extensively used in the case of a linear model. In the general case (i.e. arbitrary polynomials) a classification of the stationary points can be obtained from the linearization of the polynomial model [START_REF] Nemytskii | Qualitative Theory of Di[END_REF] but this assumes a location of critical points (i.e. to solve an arbitrary polynomial system). Instead we use the Index of a vector field to locate and characterize the stationary points. Let g = (P, Q) be a vector field defined over a Jordan curve J in the Euclidean plane, with no critical point on J. The index of g over J is proportional to the angular variation of the vector g(M) (applied at M E J) as M describes J. For the system (8), the index over an oriented Jordan curve J is given by:

Index(J) = ~ d arctan --~ p2 T Q2 (13)
A classification of the flow field g can be obtained by computing the index over a small circle surrounding an isolated critical point. Since the computation of the stationary points of system (8) leads to the problem of solving a system of polynomials with arbitrary degree, we choose to compute locally the index of g = (P, Q) over the whole flow field: At each point we consider a circle contained in the centered window W and we use the following classification:

-The index of a focus, a center or a node is equal to +l, -The index of a saddle point is -1.

Although this characterization is compendious, it characterizes the most important structures in a fluid flow field: the stationary points. The index measure, computed over all the flow, allows to obtain the critical points locations without computing the roots of the system (8). Once we locate these points we may use the linearization technique to obtain a complete description of the flow field in the neighborhood of its stationary points. Figure 3 illustrates the complete framework for processing an infrared image sequence (Meteosat) of Europe (Courtesy of LMD). A frame of the sequence is

C o n c l u s i o n

This research was done within the applicative context of environmental dynamic satellite images. We tried to solve some problem arising in processing a sequence of images representing the evolution of a physical phenomenon. This computer vision approach represents an alternative to the complex modeling of the underlying physical processes. We proposed a two stages framework allowing the processing of environmental image sequences. The first stage concerns the efficient computation of an optical flow field that preserves flow discontinuities. The later one is concerned with the interpretation of the obtained displacement field.

We are currently studying the comparison between the computer vision approach and the classical method used by oceanographic and atmospheric researchers which deal with more elaborated models.

R e f e r e n c e s

  norm of the component ~ of the motion field ~ = (u, v) along the direction of V-~. If the flow constraint equation is satisfied (i.e. dI = O) and (IVIll r O, we obtain:

  Fig.1. Example of the use of an adaptive mesh to increase the numerical accuracy of the computed flow field while reducing the algorithmic complexity of the method. We display a frame of the Hamburg Taxi sequence, the associated mesh and the optical flow field.

  ,d \ x,d \ ",,d'-,d\ "~

Fig. 3 .

 3 Fig.3. A frame of the infrared image sequence, and the plot of the computed optical flow characterizing clouds motion. The black quadrangle represents the vortex detected with the index approach. In this case, the given flow field (i.e. the optical flow) was approximated with a Q2(/f~ 2) polynomial.

  1. J.L. Barron, D.J. Fleet, and S.S. Beauchemin. Performance of optical flow techniques. International Journal of Computer Vision, 12(1):43-77, February 1994.

where P(x, y) and Q(x, y) are continuously differentiable functions. The particle trajectories are defined by the curves (~o(t), r t > 0, satisfying: (~l(t), el(t)) = (P(~(t), r

Q(~o(t), r

Modeling the orientation flow field by a dynamic system allows to characterize the flow field through the particles trajectories and their stationary points.

Different works were led on linear phase portrait models and their use for characterizing oriented texture fields [START_REF] Rao | Computerized flow field analysis: Oriented textures fields[END_REF][START_REF] Shu | Vector field analysis for oriented patterns[END_REF][START_REF] Zhong | Salient structure analysis of fluid flow[END_REF]. The main drawback of a linear phase portrait is that it can handle only one critical point. Ford and Strickland [START_REF] Ford | Nonlinear phase portrait models for oriented textures[END_REF] proposed a nonlinear phase portrait model allowing multiple critical point, but this model is computationally expensive and can not be generalized to arbitrary polynomials. This paper propose a new approach for approximating an orientation field and characterizing the stationary points of the trajectories obtained from an arbitrary polynomial phase portrait. Furthermore this model is always linear, independently of the polynomial representation.

Given two vectors A and B, a distance between them may be given by [START_REF] Golub | Matrix Computations[END_REF]:

where 0 is the angle subtended by the oriented segments. This distance represents also the area of the triangle formed by these oriented segments. Several authors [START_REF] Ford | Image models for 2-D flow visualization and compression[END_REF][START_REF] Rao | Computerized flow field analysis: Oriented textures fields[END_REF] used this measure to recover the six parameters of a linear phase model by minimizing locally the functional:

where W is the local region over-which the linear approximation is searched, 01

represents the linear model and 02 the vector field to approximate. The use the of Levenberg-Marquardt method for minimizing such a non quadratic functional, leads to a slow convergence rate [START_REF] Rao | Computerized flow field analysis: Oriented textures fields[END_REF]. In the following, we make use of another definition of the area subtended by two oriented segments that leads to the minimization of a quadratic functional. This other definition of the distance is obtained from the cross product of the vectors A and B:

Hence, considering a polynomial phase portrait model g(x,y) defined by Eq. ( 8), where P and Q E Q~(~2) = {p, st p(x,y)= Y~i,j<_,~ aljxlYJ}, we fit the model to the given orientation field ft = (fl, f2) (obtained from the optical flow field framework) by minimizing locally:

]fxg [2, (12) i,jEW

where W is a neighborhood of the image point (i,j).


Bertrand Leroy
email: bertrand.leroy@inria.fr

Isabelle L Herlin
email: isabelle.herlin@inria.fr

Laurent D Cohen
email: laurent.cohen@dauphine.fr

Face Identification by Deformation Measure

This paper studies the problem of face identification for the particular application of an automatic cash machine withdrawal: the problem is to decide if a person identifying himself by a secret code is the same person registered in the data base. The identification process consists of three main stages. The localization of salient features is obtained by using morphological operators and spatio-temporal information. The location of these features are used to achieve a normalization of the face image with regard to the corresponding face in the data base. Facial features, such as eyes, mouth and nose, are extracted by an active contour model which is able to incorporate information about the global shape of each object. Finally the identification is achieved by face warping including a deformation measure.

Introduction

Research studies related to face recognition from 2-D frontal views may be shared in two main classes according to whether they refer to a feature based approach or to a pixel approach. The first class consists in extracting the position of several facial features in order to restrict the data space. Several statistical methods, generally classical algorithms of pattern recognition, are then used to identify faces using these measures [START_REF] Yuille | Facial feature extraction by deformable templates[END_REF]. The pixel approach is a global one that requires few knowledge about face geometry. Some studies are linked to pattern matching: the measures used to identify faces are correlations between template regions. Other methods use principal component analysis applied on pixel's grey level values [START_REF] Turk | Interactive-time vision: face recognition as a visual behavior[END_REF].

The application, described in this paper, is a system similar to an automatic cash machine withdrawal where the user identifies himself by using a credit card with a secret code. It is possible to enhance the security of such a transaction by verifying that the observed face is identical to a stored face of the same person or to compare the face to information (points, curves, quantitative values, ) obtained by learning and stored in the data base.

We propose an identity check system using warping between faces as a measure of likeness. Knowing the assumed identity of a face (supposed by the secret code), the system compares this face with that one in the database corresponding to the same person. This is achieved by warping one face in order to fit the other one.

Since we have chosen to use the geometry and the locations of facial features (such as the mouth, nose and eyes) to achieve warping, each person of the database is represented by this information.

The face identification scheme consists of three main stages. In the first step, in section 2, facial features are located and the image is normalized so that features location match those of the face corresponding to this person in the database. In section 3, the features boundaries are extracted using a specific active contour model. This model is based on Fourier descriptors and is able to incorporate information about the global shape of each object. The last stage is described in section 4 and consists in measuring the warping between the normalized face and the corresponding face in the database.

Face normalization

In the context of our application, the image acquisition is achieved with a fixed lightning and is acquired from a near frontal view throughout the temporal sequence. The normalization process is based on this constraint.

The first step of our facial features localization process concerns the eyes. The nose and the mouth will be further detected using some facial geometry hypothesis.

The presence of specular spots (usually in the neighborhood of the irises) due to the direct lightning allows to locate the eyes using a "peak" morphological operator which has the effect of bringing out the local maxima of the grey level values function. A rough estimation of the eyes position is consequently obtained by looking for the two regions corresponding to the highest values of the "peak" image. Since the lightning used for the acquisition is frontal, these estimations correspond approximately to the irises positions.

In order to locate the nose and the mouth, we then use a spatio-temporal information by computing the sum of the spatio-temporal gradient norm over the sequence:

1 2 2 2 ,
where represents the grey level value, and the spatial components, the temporal component and the number of images of the sequence. is a normalization term depending on the temporal sampling of the sequence and on the motion of the different features during the sequence. As it is shown on figure 1, the image permits to pick up regions corresponding to a weak contrast but having a high deformation, like the mouth, and regions having a high contrast like the shade due to the nose relief. The localization of the nose and the mouth in the image consists in searching, along the mediatrix of the segment joining the two irises, the regions having the two higher values of the image . In order to minimize the computing time, the image is computed in small windows on both sides of the mediatrix between the two irises. Figure 2 shows the results obtained by this algorithm and it can be seen that the localization of the lower part of the nose is obtained with a good accuracy. The point corresponding to the position of the mouth depends on its motion during the sequence and is not always the same physical point. When the motion is important, this point is located within the mouth and otherwise it is generally located on the boundary of the lips.

In order to compute the image , we use the first 30 images of the sequence (with a sampling rate of 25 hz). Using a shorter sequence leads to less accurate results for the temporal information will be less significant.

In next section we make use for the system of an active contour model to extract facial features boundaries: in order to use the same initialization and parameters for the active This is achieve by applying a similarity using the irises location as references. Since the active contour parameters and the curves used for initialization are not the same for different faces, these features will be stored in the database with the facial characteristics.

Facial feature extraction

An active contour model is an energy minimizing curve guided by internal constraints and influenced by image forces that pull it towards image edges. Since the global regularity constraints used by generic active contour models (such as "snakes" [START_REF] Kass | Snakes: active contour models[END_REF][3]) are included within the internal energy function, these models do not allow to refine the space of admissible curves to a specific shape class. In order to take into account the shape of facial features, several methods based on deformable templates have been developed.

Yuille, Cohen and Hallinan [START_REF] Yuille | Facial feature extraction by deformable templates[END_REF] proposed a set of deformable templates for this purpose. Each one of these templates is designed to extract a specific feature. For example the template used to extract the eye boundary consists of a circle bounded by two parabolas and possesses 11 parameters.

Craw, Tock and Bennett [START_REF] Craw | Finding face features[END_REF] described a model which consists of a set of 40 landmark points, each one representing the position of a particular part of the structure to be located. Several algorithms are designed to deal with a particular feature using a statistical knowledge of the landmark positions. A global algorithm controls the interaction between them. A similar method has been proposed by Cootes, Tailor and Lanitis [START_REF] Cootes | Multi-resolution with active shape models[END_REF] using statistically based models which iteratively move towards structures in images after a learning process on similar images.

In this section we propose a parametrical active contour model using Fourier descriptors. The aim is to design a model sufficiently general to extract the boundary of several facial features while taking into account the specific geometry of each feature. This model needs to tune few parameters.

Actually, this model is able to extract facial features subjected to high deformation (in order to take lips motion into account) and having corners (like the mouth and the eyes).

Active contour model using Fourier descriptors

The use of Fourier descriptors for active contour models has been introduced by Staib and Duncan [START_REF] Staib | Boundary finding with parametrically deformable models[END_REF] in order to extract an object boundary. Their method is based on the use of probability distributions on the parameters. In order to be less sensitive to the initial parameter value, we propose a variational approach similar to the method used in the snake model [START_REF] Kass | Snakes: active contour models[END_REF].

An elliptic Fourier representation of a closed curve is a parametrical curve defined by:

0 cos sin ( 1
)
where is a 2 2 matrix, the number of harmonics used to describe the curve and the angular parameterization index.

Such a class of curves does not permit to represent effectively the boundaries of objects having an irregularity point within a regular section, such as the corner of the mouth or eyes. Thus, the extraction of this kind of boundary requires the use of high frequency harmonics which have the effect of creating oscillations.

where and are the parameters model, is the number of harmonics used to describe the curve and is the rotation matrix of angle . Thus, the phase shift between the two bases used for and permits to describe a half-ellipse with four parameters.

The curve modeling the boundary of the object, is obtained by minimizing a functional similar to the snake energy:

0 2 2 2
The first term inside each one of the integrals is an image potential equal to the opposite of the square of the image gradient ( 2) and the second one is an elasticity term associated to the curve tension.

In order to extract boundaries of the mouth and eyes we use one harmonic to describe the lower boundary, and harmonics for the upper boundary, . Thus has the form defined in equation ( 2) and is defined by:

Results

The algorithm used to extract features is a fine to coarse scheme (see figure 3). Defining two half-ellipses as initialization, the convergence is first achieved using a single harmonic for the curve , then the first result is used as an initialization to get a curve with two harmonics until we obtain a curve described by harmonics fitting the object boundary.

As the number of harmonics increases, the regularity constraints are smaller and the boundary can be described more precisely with a better stability.

Face Warping for identification

Several methods have been proposed to achieve face recognition using features location. Some studies use this information in order to calculate the correlation between the regions corresponding to a same feature in two different images [START_REF] Beymer | Face recognition under varying pose[END_REF]. The measure is then equal to the mean of the correlation values. Other methods utilize a more geometrical approach by defining a discriminating measure using physiological knowledge [START_REF] Kanade | Picture processing system by computer complex and recognition of human faces[END_REF], results of principal component analysis [START_REF] Costen | Automatic face recognition: what representation[END_REF] or topological graphs [START_REF] Manjunath | A feature based approach to face recognition[END_REF].

In this paper, two faces are compared by a deformation measure. Given a set of landmark points on the face, the measure is a quantification of the bending that has to be applied to the face, so that its landmark points match the ones of the other face. To thisissue, we utilize a class of deformations similar to the one used for normalization by Reisfeld [START_REF] Reisfeld | Generalized symmetry transforms: attentional mechanisms and face recognition[END_REF]. While his aim was to normalize faces in order to apply a classification on the image data, we propose to achieve the identification using the deformation measure value as the discriminating information.

Face characterization and database

The landmark points used to represent the face are located on the the features boundaries which have been extracted as described in section 3. Using all the points of the boundaries for the matching process will be useless since the information is redundant (like neighbor points on the same boundary) and some points cannot be easily extracted from the curves describing the boundaries (due to the parameterization which is not intrinsic). We made the choice of 15 points displayed on figure 4. These points characterize the position of the features as well as their global geometries.

Figure 4. points characterizing a face

However some of the chosen landmark points are mobile since they are located on a feature which is subject to temporal deformation. Therefore, in order to construct the database we use a sequence of images of people uttering a particular sentence while blinking. Then, on each image of the sequence, the features boundaries are extracted on a normalized image (as it is shown in sections 2 and 3) and for each landmark point we calculate the distribution of its position. Normality tests showed that these distributions may be represented by Gaussian laws. Thus, a person is represented in the database by a set of 4 values: where is the mean position of the landmark point and its standard deviation on the horizontal and vertical axes.

Deformation measure

Bookstein [START_REF] Bookstein | Principal warps: thin-plate splines and decomposition of deformations[END_REF] showed the interest of the thin plate splines for image registration purpose since they can be used to solve the problem of scattered data interpolation while minimizing the warping.

The total amount of bending of a surface ( being a twicely differentiable function : 2 ), can be quantified by the value of the functional: (3)

Given a set of landmark points on one face and on a second one, the function : 2 2 , that matches each point to with respect to its standard deviation and that is least bent according to the functional (3), minimizes the functional:

1 2 2 2 2
is a positive weighting parameter and is defined by: if 1 and 1 otherwise Thus, a weight is associated to each landmark point which decreases the influence of a point if its standard deviation is large. The function is constructed so that is always defined and the points having small standard deviation are not too influent.

The transformation , that minimizes , has the following form: Several types of function can be used in order to obtain a minimum for . In particular the use of 2 log (with U(0)=0) leads to an unique solution verifying the property of the thin plate spline. The coefficients 1 2 3 1 2 3 1 and 1 can be calculated by solving a linear system as described in [START_REF] Bookstein | Principal warps: thin-plate splines and decomposition of deformations[END_REF] [START_REF] Reisfeld | Generalized symmetry transforms: attentional mechanisms and face recognition[END_REF].

The quantity may be used to measure the likeness between two faces. The assumed identity will then be confirmed if this value is smaller than a given threshold. The threshold value is linked to the accuracy of the system that is wished for. If it is small, it can hapen that the identity of a person will be wrongly rejected; if it is large, the identification test will then be less restrictive. Lastely, the choice of the treshold value must take several parameters into account such as the image acquisition and the accuracy of the localization and extraction processes.

Conclusions

This paper presents a set of methods that are used to construct an identity check system. The localization of facial features (such as the eyes, nose and mouth) is achieved using a morphological operator as well as the spatio-temporal information. This method is robust provided that lightning is frontal. It is used to normalize the face image with regard to the corresponding face in the data base. The active contour model proposed for feature boundary extraction is a parametrically deformable model using Fourier descriptors. This model is designed so that few parameters are needed to extract irregular objects like the eyes or the mouth and the stability is increased by applying a coarse to fine scheme to the convergence process. A face is represented in the database by a set of landmark points belonging to the features boundaries as well as the standard deviations of the points calculated on a sequence of images. Finally the likeness between two faces is quantified using a deformation measure. This deformation corresponds to the bending to be applied on a face, so that its landmark points match the ones of the other face. To this issue, an enhancement of a standard method used to interpolate scattered data has been proposed which allows to associate a weight to each landmark point so that points having a small standard deviation have more effect on the bending.

Figure 1 .

 1 Figure 1. Opposite of the sum of the spatiotemporal gradient norm calculated on an image sequence.

Figure 2 .

 2 Figure 2. Results of the localization of the eyes, nose and mouth.

  In order to solve this problem we propose a new model made of two open curves joined at their endsuse orthogonal bases on the interval 0 , an open curve can be defined by:

  lower boundary of the nose, we use a single open curve described by equation (2).

Figure 3 .

 3 Figure 3. Initialization (first column) and results obtained with 1 (second column) and then 3 (third column) harmonics for the upper boundary of mouth and eyes and the lower boundary of nose.