Najah Hraiech
email: hraiech@gmail.com

Daniel Weinland
email: weinland@inrialpes.fr

Kamel Hamrouni
email: kamel.hamrouni@enit.rnu.tn

An Active Contour Model Based on Splines and Separating Forces to Detect the Left Ventricle in Scintigraphic Images

In this paper we present a novel method for the segmentation of scintigraphic images in order to detect the two ventricles of the heart. The method is based on the active contour model using a region based approach and parametric representation of the contour. The novel aspect of the method is that we have integrated separating forces between the two ventricles in the active contour model. The method has been applied on several scintigraphic images and compared to "active contour without edges" model. The results are satisfactory.

I.

INTRODUCTION

The need in the medical field of knowing how certain organs of the body work, especially to diagnose and determine certain heart diseases, leads to scintigraphy. Scintigraphic imaging techniques can give images of any organ of the body. The organ to be investigated is injected by a radioactive solution. A special camera measures the radiation emitted by a certain organ. One application of these techniques is to capture the geometry of the two heart ventricles.

Usually scintigraphic images are blurred with noise due to artifacts and noise from the gamma camera. As the intensity contrast between the two ventricles and the myocardium is low, the boundaries of the ventricles are often missing in certain locations. Even for skilled personal the exact boundary identification is difficult. Hence segmentation techniques can help us to recognize the two ventricles.

The aim of segmentation is to distinguish objects from background to achieve homogeneous regions. It is widely used for medical image classification.

In this paper we will present a new approach to segment the ventricles based on active contours. While existing active contour approaches can segment the heart in one "homogenous" region, most of them fail in determining the intermediate boundary, which separates the two ventricles.

For the evolution process of two active contours, centered on the right and left ventricle, we demonstrate how by integrating novel separating forces the segmentation can be efficiently improved leading to clearly separated ventricles.

The paper is organized as follows: In the first part, we review the simplified diffusion snake as introduced in [START_REF] Cremers | Diffusionsnakes: Combining Statistical Shape Knowledge and Image Information in a Variational Framework[END_REF]. Then we show how to integrate the separating forces. Finally we will give experimental results of our method in separating the two ventricles of the heart, and compare our approach with one of the most widely used active contour approaches, the Chan and Vese model [START_REF] Chan | Active Contours Without Edges[END_REF].

II. ACTIVE CONTOUR BASED ON SPLINES

The active contour model, or snake, is defined as an energy minimizing functional, where the snake's energy depends on its shape and location within the image. The goal of the active contour is to extract the boundaries of homogeneous regions that constitute an image.

The two most used active contours approaches are: Edge based segmentation approaches [START_REF] Kass | Snakes: Active contours models[END_REF] and region based segmentation approaches [START_REF] Zhu | Region competition: unifying snakes, region growing, and bayes/MDL for multiband image segmentation[END_REF][START_REF] Paragios | Geodesic active regions: A new paradigm to deal with frame partition problems in computer vision[END_REF][START_REF] Ronfard | Region-based strategies for active contour models[END_REF]. Edge based active contour approaches are based on gradient information to deform the curve towards the edges of objects within an image. Region based active contours approaches include terms describing the different regions by statistical spatial features, like mean and variance of color values. In our work we use a region based active contour approach, because scintigraphic images are often blurred with noise, and edge based segmentation is known for its sensitivity to noise.

A. Simplified diffusion snake

Our work is derived from the probably most widely studied region based Active Contour approach, which is that of Mumford and Shah [START_REF] Mumford | Boundary detection by minimizing functionals[END_REF]. Mumford and Shah proposed to approximate a given input image f with a piecewise smooth function by minimizing the following functional:

( ) ( ) , 2 1 2 1 , 2 2 2 - Ω Ω + ∇ + - = C C v dx u dx u f C u E λ (1)
where f is the input image, u is the segmented image, C the segmenting contour, and , v are regularization parameters.

Constraining the segmented image to constant regions, u(

x) = u i, for i R x ∈ , Ω ∈ i R
, the Munford-Shah functional [START_REF] Cremers | Diffusionsnakes: Combining Statistical Shape Knowledge and Image Information in a Variational Framework[END_REF] simplifies to the cartoon model as:

( ) ( ) + - = 1 0 2 2 . 2 1 , ds C v dx u f C u E s i R i i (2)
As representation for the boundary function C we choose a spline representation.

B. Spline representation

The two common ways to represent the contour function are parametric by spline curves, or implicit by the zero level set of a functional:

ℜ → Φ
, i.e. the intersection of Φ on the zero plane, see [START_REF] Sethian | Level Set Methods[END_REF].

In [START_REF] Chan | Active Contours Without Edges[END_REF] and [START_REF] Chan | A Level Set Algorithm for Minimizing the Mumford-Shah functional in image processing[END_REF] T.F. Chan and L.A. Vese explain the implementation of the Mumford-Shah functional using level sets as a representation for the contours.

In our work, we choose a parametric representation based on piecewise polynomial splines. This representation is especially suitable for the integration of the separating forces, which are directly attached onto the control points of the splines.

For the experiments we used the closed, uniform, quadratic and periodic B-spline curve, for more detail see [START_REF] Farin | Curves and Surfaces for Computer-Aided Geometric Design[END_REF] and [START_REF] Bartels | An introduction to Splines for use in Computer Graphics and Geometric Modeling[END_REF]. Using this spline approximation the contour C has the following form:

[ ] , 1 , 0 : Ω → C ( ) ( ), 1 s B p s C n N n n = = ( 3
)
where B n are the uniform, periodic, quadratic B-spline basis functions and p n = (x n, y n ) T denote the control points.

C. Curve Evolution

We minimize the energies of the simplified diffusion snake (2), using gradient descent with respect to both, the segmenting contour C and the segmented image u, what leads to the following Euler-Lagrange equation:

( ) ( ) [ ] ( ) ( ) 0 = ′ ′ - ⋅ - = + - s C v s n s e s e dC dE [ ], 1 , 0 ∈ ∀s (4)
where n is the outer normal vector on the contour for each control point. The terms e + and e -denote the energy density outside and inside the contour C(s) respectively, see [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF], defined by:

( ) . 2 u f e - = - + (5)
We solve the minimization problem by gradient descent. The result is summarized in equation ( 6): . s ∀ (6) In order to obtain an evolution rule for the control points of the spline C we insert the spline curve (3) in equation ( 6) to derive the following evolution equation for the control points:

( ) ( ) ( ) ( ) [ ] ( ) ( )
( ) ( ) ( ) ( ) [ ] ( ) ( ) ( ) . , , , 2
N i i i N i i = - + = + ⋅ - = (7)

D. Implementation

Discretizing equation ( 7) with a set of nodes s i along the contour gives us a set of linear differential equations. The solution gives the temporal evolution for the coordinates of each control point (x m , y m ):

( ) ( ) ( ) ( ) [ ] , 2 1 1 1 1 + - - + = - + - + ⋅ - = i i i x s s N i mi m x x x v n e e B dt t dx i i (8) ( ) ( ) ( ) ( ) [ ] . 2 1 1 1 1 + - - + = - + - + ⋅ - = i i i y s s N i mi m y y y v n e e B dt t dy i i (9)
The matrix B contains the spline basis functions evaluated at the nodes s i : B i,j =B i (s j ), ( s i corresponds to the maximum of B i ). Notice, that the matrix B is cyclic and tridiagonal.

The first term in the set of equations ( 8) and ( 9) maximizes the homogeneity in the adjoining regions as measured by the energy densities [START_REF] Paragios | Geodesic active regions: A new paradigm to deal with frame partition problems in computer vision[END_REF]. The homogeneity in the adjoining regions forces the contour to recognize the boundaries of the object. The second term smoothes the contour and enforces an equidistant spacing between successive control points.

In our implementation we need at the beginning the initialization of two set of control points (each set corresponds to one ventricle), then we interpolate them by the spline contour. After that, we implement the evolution of this curve as explained in equations ( 8) and [START_REF] Chan | A Level Set Algorithm for Minimizing the Mumford-Shah functional in image processing[END_REF].

Note, that with each new iteration, we have first to calculate the new control points (the previous points + the evolution) and then calculate new values for e + and e -depending on the new splines. These two steps are iterated until the difference between the subsequence evolutions is certainly small. For equation [START_REF] Paragios | Geodesic active regions: A new paradigm to deal with frame partition problems in computer vision[END_REF] we take the difference between inside mean and value of the pixel lying on the inside pointing unit normal of the control; and outside respectively.

E. Discussion

Applying the described actives contours to sample images we realize that the results depend on the initialization: when we have a good initialization the active contour recognizes the two ventricles. With random initialization the results were not acceptable as the two ventricles are often not separated, illustrated in figure 1 and figure 2.

In the next section we will integrate separating forces in the curve evolution in order to prevent this effect.

III. INTEGRATION OF SEPARATING FORCES BETWEEN THE TWO VENTRICLES

Usually active contours based on splines are a suitable model to recognize boundaries within an image. Nevertheless, using one active contour for every ventricle, experiments showed their tendency to merge ignoring the separating boundary. Therefore, we suggest introducing additional forces into the curve evolution in order to keep the two splines separated from each other.

These forces depend on the distance between the control points on the two contours and are defined as: , exp so that F j is the force from all control points i onto the control point j , Adding the forces ( 10) and ( 11) into equations ( 8) and ( 9) leads to the following temporal evolution equation for the coordinates of each control point (x m , y m ): Regularization parameter h>0 determines the trade-off between the segmenting and the separating forces.

( ) ( ) ( ) ( ) [ ] , 2 1 1 1 1 x i i i i x s s N i mi m n F h x x x v n e e B dt

IV. EXPERIENTAL RESULTS

To evaluate our work, we test our implementation with scintigraphic images of different quality. In our implementation we manually initialize the splines. Hereby the number of control points is of importance (see below), while small variations in the initial shape did not affect the results. For the regularization parameters we use an empirical choice.

For most of the images, the tendency of the two splines to merge (figure 1) was successfully suppressed and the evolution of the contours stabilizes around the true boundary of the ventricles (figure 3). Figure 4 shows a comparison of our results with that derived from the Chan and Vese approach [START_REF] Chan | Active Contours Without Edges[END_REF]. Our method showed to be superior in separating the two ventricles (e.g. figure 4 top row), but has limitations depending on the number of control points (figure 5).

The limitation for the control points is demonstrated as follows: We start from a small contour and the contour grows to cover the boundaries, but has not a sufficient number of control points to describe the exact shape of the ventricles (figure 5).

For further work we suggest to use a dynamical number of control points, to circumvent the later problem. V. CONCLUSION

In this paper, we proposed a new active contour method using a uniform quadratic B-spline linked with extra forces to implement a region-based active contour segmentation algorithm for scintigraphic images.

We conclude that the integration of separating forces into active contours model improves the segmentation of the two ventricles, compared to common active contour approaches.

To show the performance of our results, we compared it with the "Active Contour without edges" approach [START_REF] Chan | Active Contours Without Edges[END_REF].

Figure 1 .

 1 Figure 1. Result of segmentation using Active Contour based on Splines. False initialisation leads to poor final results: (left) original image, (middle) evolution of spline curves, (right) the two curves merge before the splines find the true boundaries.

Figure 2 .

 2 Figure 2. Result of segmentation using Active Contours based on Splines using good initialization: (left) original image, (middle) evolution of spline curves, (right) the segmented images.

  so that F i is the force from all control points j onto the control point i ,

  The variance σ controls the region of influence of the forces. distance between two controls points belonging to different curves.

Figure 3 .

 3 Figure 3. Results using integration of Active Contours based on Splines and forces between the regions of interest: (left) the initialization of spline, (middle) the evolution of spline curves, (right) the segmented images.

Figure 4 .Figure 5 .

 45 Figure 4. Results of segmentation of the image of the heart: (left) original images, (middle) segmentation with "Integration of separating forces to Active Contours based on Splines", (right) segmentation of original images with "Active Contours Without Edges". (top) the Chan and vese approach fails to separate the two ventrices. (bottom) variations in segmentation dependent on contpour representation (spline vs. level-set)

ACKNOWLEDGMENT

This work was partially done while visiting the CVGPR group, University of Mannheim, Germany.


Fabio Cuzzolin

Stefano Soatto

Learning Riemannian Metrics for Classification of Dynamical Models

Consider the problem of classifying motions, encoded as dynamical models of a certain class. Standard nearest neighbor classification then reduces to find a suitable distance function in the space of the models. In this paper we present a supervised differential-geometric method to learn a Riemannian metric for a given class of dynamical models in order to improve classification performances. Given a training set of models the optimal metric is selected among a family of pullback metrics induced by the Fisher information tensor through a parameterized diffeomorphism. Experimental results concerning action and identity recognition based on simple scalar features are shown, proving how learning a metric actually improves classification rates when compared with Fisher geodesic distance and other classical distance functions.

Introduction

Human motion recognition is one of the most popular fields in computer vision, due to both its applicative potential and its richness in terms of the technical issues involved. Consider then the problem of classifying a number of movements, represented as sequences of image features: Representing those sequences in a compact way as simple dynamical models has proved to be effective in problems like dynamic textures or gait identification. The motion classification problem then reduces to finding the appropriate distance function in the space of the dynamical models of the chosen class. A number of distance functions for linear systems has been already introduced, in particular in the context of system identification: Martin's distance between cepstrums [START_REF] Martin | A metric for ARMA processes[END_REF], subspace angles [START_REF] Cock | Subspace angles and distances between arma models[END_REF], gap metric [START_REF] Zames | Unstable systems and feedback: The gap metric[END_REF] and its variants nu-gap [START_REF] Vinnicombe | A ν-gap distance for uncertain and nonlinear systems[END_REF] and graph metric [START_REF] Vidyasagar | The graph metric for unstable plants and robustness estimates for feedback stability[END_REF], kernel methods [START_REF] Smola | Hilbert space embeddings in dynamical systems[END_REF]. Besides, a vast literature can be found about dissimilarity measures between hidden Markov models [START_REF] Silva | Average divergence distance as a statistical discrimination measure for hidden Markov models[END_REF][START_REF] Lyngso | Metrics and similarity measures for hidden Markov models[END_REF], most of them variants [START_REF] Do | Fast approximation of Kullback -Leibler distance for dependence trees and hidden Markov models[END_REF] of the Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF][START_REF] Juang | A probabilistic distance measure for hidden Markov models[END_REF]. However, a simple mental experiment is enough to understand how no single distance function can possibly outperform the others in each and every classification problem, since labels can be assigned arbitrarily to any given dataset.

The most reasonable thing to do, when possessing some a-priori information in terms of partially labelled data or similarity classes, is then try and learn in a supervised fashion the "best" distance function for a precise classification problem. This topic has become quite popular in the last few years (see for instance [START_REF] Bar-Hillel | Learning distance functions using equivalence relations[END_REF][START_REF] Bar-Hillel | Learning distance functions for image retrieval[END_REF][START_REF] Bilenko | Integrating constraints and metric learning in semi-supervised clustering[END_REF][START_REF] Schultz | Learning a distance metric from relative comparisons[END_REF][START_REF] Tsang | Distance metric learning with kernels[END_REF][START_REF] Zhang | Learning metrics via discriminant kernels and multidimensional scaling: Toward expected euclidean representation[END_REF][START_REF] Eick | Using clustering to learn distance functions for supervised similarity assessment[END_REF]). Many unsupervised algorithms, in particular, take an input dataset and embed it in some other space, implicitly learning a metric (locally linear embedding [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] among the others), but fail to learn a full metric for the whole input space. In this paper we propose a differential-geometric method that, given a dataset of unlabeled linear systems of a given class, allows to learn the Riemannian metric that minimizes the inverse volume element around the data, this way forcing the geodesic of the manifold to pass through the region where the data are actually distributed. We will show how this improves the classification performance by posing the problem in the region where the systems actually live. More precisely, if the models belong to a Riemannian manifold M , any diffeomorphism of M onto itself induces a pullback metric on M . By designing a suitable family of diffeomorphisms depending on a parameter p we then obtain a family of pullback metrics on M . We can then analytically compute the volume element of the induced metric as a function of p, and pose a well-defined optimization problem to get the value of the parameter.

To prove the advantages of a classification scheme based on a metric learnt from the data, we consider two different problems: action recognition and identity recognition from gait. We use the image sequences collected in the Mobo database [START_REF] Gross | The cmu motion of body (mobo) database[END_REF] to show experiments in which simple nearest-neighbor classifiers based on the pullback of the classical Fisher metric between stochastic models outperform analogous classifiers based on all the known distances between models.

The study of the geometrical structure of the space formed by a family of probability distribution is first due to Rao [START_REF] Rao | Information and accuracy attainable in the estimation of statistical parameters[END_REF], and was developed by Nagaoka and Amari [START_REF] Amari | Differential geometry of curved exponential families: curvatures and information loss[END_REF][START_REF] Amari | Differential geometric methods in statistics[END_REF]. A family S of probability distributions p(x, ξ) depending on a n-dimensional parameter ξ can be regarded in fact as an n-dimensional manifold. If the Fisher information matrix

g ij . = E ∂logp(x, ξ) ∂ξ i ∂logp(x, ξ) ∂ξ j is nondegenerate, G = [g ij ]
is a Riemannian tensor, and S is a Riemannian manifold. The squared distance ds 2 between two nearby distributions p(x, ξ) and p(x, ξ + ∂ξ) is known to be twice the respective Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF]. Some approximations of this measure have been introduced based on Monte-Carlo methods or simplifying approximations [START_REF] Juang | A probabilistic distance measure for hidden Markov models[END_REF][START_REF] Silva | Average divergence distance as a statistical discrimination measure for hidden Markov models[END_REF].

The problem of defining a metric for linear dynamical systems has been studied also in the framework of robust control.

One criterion is to compare the outputs when the inputs are restricted to the class of inputs that give bounded outputs. This approach was introduced by Zames and El-Sakkary [START_REF] Zames | Unstable systems and feedback: The gap metric[END_REF] [24] using the notion of gap metric. The original rationale for the (linear) gap metric was to provide a suitable topology in which small errors in the gap in open-loop systems would correspond to small errors in the norm of the stable closed-loop. Another solution called ν-gap metric to the problem of comparing two systems that is appropriate for feedback analysis was suggested by Vinnicombe [START_REF] Vinnicombe | A ν-gap distance for uncertain and nonlinear systems[END_REF].

Another distance function has recently been introduced based on the notion of cepstrum, the inverse Fourier transform of the logarithm of the spectral density. Martin [START_REF] Martin | A metric for ARMA processes[END_REF] defined a new metric for the set of SISO ARMA models, based on the inner product of the cepstra

(c n ) n∈Z , c * n = c -n , yielding a Euclidean metric d(c, c ) = |c -c | 2 = n n|c n -c n | 2 . De Cock [2]
showed how the cepstrum norm of a SISO ARMA model is a function of the principal angles between the row space of the controllability matrices of the model and its inverse.

Learning metrics from the data: pullback metrics

None of the above distances (which have usually been designed for a precise goal) can obviosly outperform the others in every classification problem. Labels can be assigned arbitrarily to dynamical systems, so that a distance which works well for a particular labeling would fail if employed in a different classification task. When having some knowledge about the dataset (similarity between pairs, partial labeling, etc.) it makes much more sense to learn a distance metric reflecting this information, and use the learnt metric for classification (see [START_REF] Xing | Distance metric learning with applications to clustering with side information[END_REF]). Many unsupervised algorithms take an input dataset and embed it in some other space, implicitly learning a metric ( [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] among the others). However, they fail to learn a full metric for the whole input space, so that the embedding has usually to be recomputed when new data are available 1 .

Some simple notions of differential geometry can provide us with a fundamental tool to build a structured family of metrics on which to define an optimization problem, the basic ingredient of a metric learning algorithm. Consider a family of diffeomorphisms between the Riemannian manifold M in which the dataset D ⊂ M resides and itself:

F λ : M → M, m → F λ (m) (1)
This family of maps induces a family of pullback metrics on the space itself. Let us call T m M the tangent space to M in m. Any diffeomorphism F is associated with a push-forward map

F * : T m M → T F (m) M v ∈ T m M → v ∈ T F (m) M defined as (F * v)f = v(f • F )
for all the smooth functions f on M . Then, given a Riemannian metric g : T M × T M → R on M , the diffeomorphism F induces a pullback metric on M as follows:

g * m (u, v) . = g F (m) (F * u, F * v). ( 2
)
Any parametric family of differentiable maps (1) generates naturally a parametric family of metrics on the original space M . The geodesic (shortest path) connecting two points under the pullback metric is the lifting of the geodesic associated with the original metric. In other words, the pullback geodesic between two points is just the geodesic connecting their images with respect to the original metric.

Volume element minimization

As our task is to classify data, a natural optimization criterion would be to maximize the classification performance achieved by using the learnt metric for a simple nearest-neighbor classifier. Xing et al. [START_REF] Xing | Distance metric learning with applications to clustering with side information[END_REF] have recently proposed a way to solve this optimization problem for linear maps y = A 1/2 x, when some pairs of points are known to be "similar". This leads to a parameterized family of Mahalanobis distances x -y A . Given the linearity of the map the problem of minimizing the squared distance between similar points turns out to be convex, hence solvable with efficient numerical algorithms. Other people [START_REF] Shental | Adjustment learning and relevant component analysis[END_REF] have successfully faced the same problem in a linear framework, by solving an optimization problem based on information theory (relevant component analysis).

Unfortunately, maximizing classification performance is in practice very hard in a nonlinear context. A reasonable approach consists then on choosing a different objective function, which has to be correlated with the natural one in the sense that it should improve classification performances. An interesting choice has been recently suggested by G. Lebanon [START_REF] Lebanon | Computing the volume element of a family of metrics on the multinomial simplex[END_REF] in the context of document retrieval (see also [START_REF] Murray | Differential Geometry and Statistics[END_REF]). The idea is to minimize the volume element associated with the pullback metric, in order to force the geodesics to pass through the most populated regions of the original space. More precisely, the function to minimize is

O(D) = N k=1 (det g(m k )) -1 2 M (det g(m)) -1 2 dm . ( 3
)
where g(m k ) denotes the Riemannian metric at the point m k of a dataset living on a Riemannian manifold M . This is clearly related to finding a lower dimensional representation of the dataset, in a similar fashion to dimensionality reduction techniques like locally linear embedding [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] or laplacian eigenmaps [START_REF] Belkin | Semi-supervised learning on riemannian manifolds[END_REF] 2 . The function [START_REF] Zames | Unstable systems and feedback: The gap metric[END_REF] to optimize can also be seen as the inverse of Jeffreys' prior [START_REF] Jeffreys | Theory of Probability[END_REF][START_REF] Lebanon | Computing the volume element of a family of metrics on the multinomial simplex[END_REF].

Algorithm

To find the expression of the Gramian det g * associated with the pullback metric (2) we first need to choose a base of the tangent space T m M to M . Let us then denote with {∂ i }, i = 1, ..., dim M the base of T m M . We can then compute the push-forward of the vectors of this base, getting in turn a base for T F (m) M . By definition [START_REF] Petersen | Riemannian Geometry[END_REF], the push-forward of a vector

v ∈ T m M is F p (v) . = d dt F p (m + t • v) t=0 , v ∈ T m M. ( 4
)
The diffeomorphism F p then induces a base for the space of vector fields on M , w i . = {F p (∂ i )}, for i = 1, ..., dim M . We can rearrange these vectors as rows of a matrix J = [w 1 ; • • • ; w dim M ]. The volume element of the pullback metric g * in a point m ∈ M is then given by the determinant of the Gramian [START_REF] Lebanon | Computing the volume element of a family of metrics on the multinomial simplex[END_REF],

det g * (m) . = det[g(F * p (∂ i ), F * p (∂ j ))] ij = det(J T gJ).
If J is a square matrix (as in the following) we get simply

det g * (m) = det(J) 2 • det g(m).
We can then easily find the expression of the objective function [START_REF] Zames | Unstable systems and feedback: The gap metric[END_REF]. Of course, once found the optimal metric, standard classification algorithms require to measure the geodesic distance between any two points. As geodesics for g * p are liftings, we only need to know the geodesics of the original Riemannian space.

Riemannian manifold structure under Fisher metric [START_REF] Rijkeboer | Differential geometric models for time-varying coefficients of autoregressive processes[END_REF]. A natural parametrization uses the non-unit coefficients (a 1 , a 2 ) of the denominator of the transfer function,

h(z) = z 2 /(z 2 + a 1 z + a 2 ) (which corresponds to the AR difference y(k) = -a 1 y(k -1) -a 2 y(k -2)).
To impose stability the necessary conditions are

1 + a 1 + a 2 > 0, 1 -a 1 + a 2 > 0, 1 -a 2 > 0.
The manifold is then composed by a single connected component (see Figure 1). The Riemannian Fisher tensor can be expressed as [START_REF] Peeters | On the Riemannian manifold structure of classes of linear systems[END_REF] g

(a 1 , a 2 ) = 1 (1 + a 1 + a 2 )(1 -a 1 + a 2 )(1 -a 2 ) 1 + a 2 a 1 a 1 1 + a 2 . ( 5
)
However, alternative local coordinates are given by the Schur parameters γ 1 = a 1 /(1 + a 2 ), γ 2 = a 2 under which the metric tensor simplifies as

g(γ 1 , γ 2 ) = 1 (1 -γ 2 2 ) (1+γ 2 ) 2 (1-γ 2 1 ) 0 0 1 (6)
while the parameter domain becomes Consider now the class of stable discrete-time linear SISO systems of order 1,

|γ 1 | < 1, |γ 2 | < 1. [0,-1] [2,1] [-2,1]
x(k + 1) = ax(k) + bu(k) y(k) = cx(k) ( 7
g(r, θ) = 1 0 0 r 2 (8)
and the manifold M (1, 1, 1) coincides with the two-dimensional plane, vertical axis excluded (r = 0).

Geodesics

We can use the analytic expressions of the Fisher metric to write the geodesic equations for these classes of linear systems. Using Einstein's notation the geodesic equation reads

d 2 dt 2 x a + Γ a bc d dt x b d dt x c = 0
where {Γ a bc } are the Christoffel coefficients of the second kind [START_REF] Petersen | Riemannian Geometry[END_REF], Γ a bc = Γ bcd g ad where g ad = g -1 (a, d) are the entries of the inverse of the metric tensor, and the Christoffel coefficients of the first kind are a function of the derivatives of the metric:

Γ abc . = 1 2 g ad ( d dx c g bd + d dx b g cd - d dx d g cb ).
Ravishanker [START_REF] Ravishanker | Differential geometry of ARMA models[END_REF] has explicitly computed the expression for the Fisher metric of ARMA(p,q) systems using as parametrization the inverses of poles and zeros, and we could then solve the geodesic differential equation analytically. However, all the geodesics of stable AR(2) systems endowed with the Fisher metric (6) as a function of the Schur parameters have been analytically computed by Rijkeboer [START_REF] Rijkeboer | Fisher optimal approximation of an AR(n)-process by an AR(n-1)-process[END_REF][START_REF] Rijkeboer | Differential geometric models for time-varying coefficients of autoregressive processes[END_REF]. The analytical expressions of the geodesics of the manifold M (1, 1, 1) endowed with the Riemannian tensor induced by the scalar product of a Hilbert space H 2 have been instead computed by Hanzon [START_REF] Hanzon | Identifiability, recursive identification and spaces of linear dynamical systems[END_REF] through an isometry with a Riemannian surface. We will use it instead of the Fisher tensor in the following 3 . More precisely, given two systems Σ 1 = (r 1 , θ 2 ) and Σ 2 = (r 2 , θ 2 ) their geodesic distance according to the metric (8) is

d(Σ 1 , Σ 2 ) = (r 2 1 -2r 1 r 2 cos(θ 2 -θ 1 ) + r 2 1 ) 1/2 if |θ 1 -θ 2 | ≤ π and sgn(r 1 ) = sgn(r 2 ); d(Σ 1 , Σ 2 ) = |r 1 | + |r 2 |
otherwise.

Volume elements of pullback metrics for spaces of linear models

The algorithm of Section 3.1 adapts easily to the case in which the dataset is a collection of linear systems.

One natural choice for a diffeomorphism of AR( 2) onto itself is suggested by the simplicial form of the manifold (see Figure 1),

F p (m) = F p ([m 1 , m 2 , m 3 ]) = 1 λ • m [λ 1 • m 1 , λ 2 • m 2 , λ 3 • m 3 ] (9)
where

p = λ = [λ 1 , λ 2 , λ 3 ] with λ 1 + λ 2 + λ 3 = 1, while m = [m 1 , m 2 , m 3 ]
collects the simplicial components of a system a in the manifold:

a = [a 1 , a 2 ] = m 1 [0, -1] + m 2 [2, 1] + m 3 [-2, 1] , a ∈ AR(2). ( 10
)
Theorem 1 The volume element of the pullback metric on AR(2) induced by the diffeomorphism ( 9) is given by

detg * p (m) ∝ (λ 1 λ 2 λ 3 ) 2 (λ • m) 6 • 1 m 2 1 m 2 m 3 .
Let us choose as base of the tangent space in AR(2) the set

∂ 1 = [1/2, 1/2] , ∂ 2 = [-1/2, 1/2] .
To compute the Gramian we need to express the diffeomorphism with respect to a. From Equation ( 10)

a 1 = 2(m 2 -m 3 ), a 2 = m 2 + m 3 -m 1 , while the inverse relation is m 2 = 1+a 1 +a 2 4 , m 3 = 1-a 1 +a 2 4 , m 1 = 1-a 2 2 . Hence F λ can be expressed in terms of a 1 , a 2 as F λ (a) = 1 ∆ [2λ 2 (1 + a 1 + a 2 ) -2λ 3 (1 -a 1 + a 2 ), λ 2 (1 + a 1 + a 2 ) + λ 3 (1 -a 1 + a 2 ) -2λ 1 (1 -a 2 )]
where

∆ = 2λ 1 (1 -a 2 ) + λ 2 (1 + a 1 + a 2 ) + λ 3 (1 -a 1 + a 2 ), so that F λ (a + tv) = [2λ 2 (1 + a 1 + tv 1 + a 2 + tv 2 ) -2λ 3 (1 -a 1 -tv 1 + a 2 + tv 2 ), λ 2 (1 + a 1 + tv 1 + a 2 + tv 2 ) + λ 3 (1 -a 1 -tv 1 + a 2 + tv 2 ) -2λ 1 (1 -a 2 -tv 2 )] .
We can then compute4 d dt F λ (m + t • v) t=0 , and in particular

w 1 = d dt F λ (m + t • ∂ 1 ) t=0 = [2λ 1 λ 2 (3 -a 2 + a 1 ) + 2λ 3 (2λ 2 -λ 1 )(1 -a 1 + a 2 ), 2λ 1 λ 2 (3 -a 2 + a 1 ) + 2λ 1 λ 3 (1 -a 1 + a 2 )]
while

w 2 = d dt F λ (m + t • ∂ 2 ) t=0 = [-2λ 1 λ 3 (3 -a 2 + a 1 ) + 2λ 2 (λ 1 -2λ 3 )(1 + a 1 + a 2 ), 2λ 1 λ 3 (3 -a 2 + a 1 ) + 2λ 1 λ 2 (1 + a 1 + a 2 )].
The determinant of the matrix J is then (after a few passages)

det J = 32 λ 1 λ 2 λ 3 ∆ 3 = 1 2 λ 1 λ 2 λ 3 (λ • m) 3 .
Eventually, the function (3) to maximize assumes the form

O(p) = N k=1 (λ • m k ) 3 AR(2) (λ • m) 3 m 1 √ m 2 m 3 dm (11)
where the normalization factor

I(λ) = AR(2) (λ • m) 3 m 1 √ m 2 m 3 dm
forbids trivial solutions in which the volume element is minimized by shrinking the whole space. It can be computed by decomposing the cube (λ • m) 3 using Tartaglia's formula, obtaining

I(λ) = c 1 +c 2 +c 3 =3 3! c 1 !c 2 !c 3 ! 3 j=1 λ cj j AR(2) m 1+c 1 1 m 1/2+c2 2 m 1/2+a3 3 dm.
About M (1, 1, 1), Henon's map is the simplest non-trivial diffeomorphism of the plane: x = y + x 2 + a, y = -bx. Here we choose a modified version:

F p (r, θ) = [ bθ + ar 2 + br, br + aθ 2 + bθ ], (12)
with p = [a, b].

Theorem 2

The volume element of the pullback metric on M (1, 1, 1) induced by the diffeomorphism ( 12) is given by

detg * p (x) = 4a 2 (2arθ + br + bθ) 2 r 2 .
It is easy to realize that in this case

J = 2ar + b b b 2aθ + b
so that det J = 2a(2arθ + br + bθ), while from Equation (8) det g(r, θ) = r 2 . This time the function (3) can be expressed as

O(D) ∝ k 1 r k (2ar k θ k ) + br k + bθ k • 1 I(a, b) ( 13
)
where I(a, b) is again the normalization factor. Note that when I(p) is hard to compute, or the original manifold has no finite volume (and hence does not admit Jeffrey's prior), an acceptable replacement is

I(p) = N k=1 (detg(x)) -1/2
. The functions ( 11) and ( 7) can be maximized by means of any numerical optimization scheme, simple (like gradient descent) or more sophisticated.

Human motion experiments

To test the effectiveness of the approach we considered a significant vision application, namely the recognition of actions and identities from image sequences. We used the Mobo database [START_REF] Gross | The cmu motion of body (mobo) database[END_REF], a collection of 600 image sequences of 25 people walking on a treadmill in four different variants (slow walk, fast walk, walk on a slope, walk carrying a ball), seen from 6 different viewpoints corresponding to cameras equally distributed around the subject (see Figure 2-left). Each sequence of the database possesses three labels: action, view, and identity.

The database comes already with preprocessed silhouettes of the moving body as black and white images. As the results of the above sections concern only single input-single output linear systems, we chose to extract from each image the width of the minimum box containing the silhouette. This way, each image sequence is represented as a scalar signal. Each scalar sequence has been then passed as input to an identification algorithm: we hence generated a dataset of linear systems of the two classes discussed above, one system for each labeled motion sequence.

This metric learning approach is based on the assumption that the two optimization problems, i.e. maximizing classification performance and minimizing volume elements are indeed correlated. We empirically evaluated the goodness of this conjecture by measuring the performance of a nearest-neighbor classifier based on the optimal pullback metrics induced by the diffeomorphisms ( 12), ( 9), and comparing the results with those of NN classifiers based on all the other known distances (included the Fisher geodesic distance).

Identity recognition

In a first experiment we selected a training set of models, and used the pullback geodesic distance to classify the identity of the person appearing in a different set of randomly selected sequences. This is a very difficult problem, as there are 25 different people, and the one-dimensional signal we chose to represent sequences clearly provides insufficient information. However, measuring the comparative performance of the metrics can be useful to see how learning a metric appositively gives an advantages when using a NN classifier. We implemented a naive Frobenius norm of the system matrices in canonical form, gap and ν-gap metrics [START_REF] Zames | Unstable systems and feedback: The gap metric[END_REF][START_REF] Vinnicombe | A ν-gap distance for uncertain and nonlinear systems[END_REF], subspace angles [START_REF] Cock | Subspace angles and distances between arma models[END_REF][START_REF] Martin | A metric for ARMA processes[END_REF], and of course the Fisher geodesic distance together with the associated optimal pullback metric.

Figure 2-right shows the percentage of correctly classified testing systems over several runs in which we randomly selected an increasing number of systems in both training and testing set. We also tested a different diffeomorphism for AR(2) systems, namely

F λ (m) = [ λm 1 + (1 -λ)m 2 , λm 2 + (1 -λ)m 3 , λm 3 + (1 -λ)m 1 ], (14)
with 0 < λ < 1. We can notice how the other distances do not show any significantly different behavior. This of course was to expect, since they have not been designed to solve classification tasks (usually for identification purposes). This has been the case in all our experiments. Hence, in the following we will only report the performance of the second best distance for comparison. In a second experiment, we selected a training set of models, and used the pullback geodesic distance to classify the action performed by the person in a different set of randomly selected sequences. Recall that all the actions in the Mobo database are just slightly different variations of the walking gait, and we were using a single scalar feature. Figure 3 illustrates the relative performances of the optimal pullback metric for autoregressive systems under deformation map [START_REF] Bilenko | Integrating constraints and metric learning in semi-supervised clustering[END_REF] and its best competitor. The correct classification rate has been measured for different sizes of both training and testing set, by randomly selecting for each class of action an equal increasing number of systems to build the two collections. To test the approach more thoroughly we also conducted six separate experiments by selecting the portion of the dataset associated with a single view, for each possible view. These are completely different classification problems, as the datasets can turn out to live in different submanifolds. It is in fact interesting to try and see what happens when we choose a different class of diffeomorphisms. Figure 4 illustrates the average performance of the classifiers associated with the optimal pullback metric for AR(2) models with (this time) diffeomorphism ( 9) and the best competing distance for four view-dependent experiments. This time increasing abscissae (from 1 to 7) mean an increased testing set size. As usual the average is computed by repeated random selections, and the optimal metric performs far better than the others.

Of course the choice of the family of deformations affects the range of candidate metric tensors the algorithm can select from. We have seen as a matter of fact how the first family of mappings, which depends on 2 parameters instead of 1, produces a much wider selection choice for the optimization algorithm, eventually delivering a better performance. The same phenomenon can be appreciated in the identity experiment (see Figure 2). It is really interesting to understand whether the size of the training set we use to learn the optimal metric has or not an influence on the classification performance. Of course, as our method exploits the training set to understand the local structure of the data, it is reasonable to conjecture that larger training sets should lead to better recognition rates, as long as the training data better represent the whole dataset. This is confirmed by Figure 5, where the recognition rate is plotted against the size of the training set. Finally, let us visualize the effect of the optimal diffeomorphism on the dataset. We can easily do it in the action recognition experiment, as there are only 4 labels. Figure 6-left shows an instance of the original dataset represented as a set of points in the AR(2) simplex. Points of the training set are drawn as squares, while systems in the testing set are represented as crosses. Figure 6-right shows how the diffeomorphism deforms the space, plotting the images of the dataset according to the optimal map. Data-points are colored according to their true class. We can visually appreciate how the deformation improves the separation between clusters.

Conclusions

In this paper we discussed the problem of learning the "best" metric for a classification problem involving linear dynamical models. Given a training set of models we pose a related optimization problem in which the pullback metric induced by a diffeomorphism which minimizes the volume element around the available data is learned. We adopt as basis metric tensor the classical Fisher information matrix. This yields a global embedding, while usual spectral methods only provide images of training points.

We have also shown experimental results concerning identity and action recognition, which proved how such a learnt metric actually improves classification performance with respect to the Fisher metric and all the other known distances between dynamical models. Setting and solving the minimization problem for MIMO systems is the next natural step. This will be of course of paramount importance to test this approach in a more realistic scenario.

Figure 1 :

 1 Figure 1: Graphical representation of the manifold of stable autoregressive systems of order 2, AR(2), with the non-unit coefficients of the denominator of h(z) as parameters. It forms a simplex with vertices [-2, 1], [2, 1], [0, -1].

√ 1 -a 2 ,

 12 ) for which we can choose a canonical representation by setting c = 1, so that the transfer function is h(z) = b/(z -a). We denote this class of systems with M (1, 1, 1). Under the conditions |a| < 1 (stability) and b = 0 (minimality) this family of linear models form a manifold with local coordinates (a, b), consisting of the two connected components {(a, b) : a < 0, b < 0} and {(a, b) : a < 0, b > 0}. Using the generalized polar coordinates r . = b θ =atanh(a) the Fisher tensor becomes

Figure 2 :

 2 Figure 2: Left: location and orientation of the six cameras in the Mobo experiment (the origin of the frame is roughly in the position of the walking person on the treadmill). Right: Performance of NN classifier based on all the known distances between dynamical models, included the pullback metric induced by the Fisher geometry. Correct classification rate for the used metrics averaged on the size of the training set, in the ID experiment. 1 -Frobenius norm of the matrices in canonical form; 2 -gap metric; 3-ν-gap matric; 4 -subspace angles; 5 -basis Fisher metric; 6 -pullback metric with diffeomorphism (14); 7 -pullback metric with diffeomorphism (9).

Figure 3 :

 3 Figure 3: Top: Recognition performance of second-best distance (left) and optimal pullback metric for AR(2) systems and diffeomorphism (14) (right) in the action recognition experiment. Here the whole dataset was considered, regardless of the viewpoint. Increasing x indices (1 : 50) mean decreasing size of the testing set. Increasing y indices (1 : 5) represent increasing size of the training set. Bottom-left: Difference between the two performances. Bottom-right: relative performances for view 5 only. Again, the pullback metric outperforms the best competitor.

Figure 4 :

 4 Figure 4: Correct action classification rates for four different viewpoints. Sequences are represented as autoregressive order 2 systems. The performance (in red) of the pullback metric associated with the map (9) is shown outperforming the best competitor (in blue) for increasing sizes of the testing set.

Figure 5 :

 5 Figure 5: Effect of the size of the training set on the recognition accuracy in the action recognition experiment, M (1, 1, 1) systems. As the size of D increases (abscissa indices) the number of correct matches grows in both the overall experiment for all view (left) and for view 2 (right).

Figure 6 :

 6 Figure 6: Effect of the optimal diffeomorphism (9) on the dataset.

Even though some work on out-of-sample extensions of spectral methods has been done[START_REF] Bengio | Out-of-sample extensions for LLE, Isomap, MDS, eigenmaps, and spectral clustering[END_REF].

Learning metrics for linear dynamical modelsEven though large classes of non-linear models can be endowed with a Fisher metric tensor, most efforts have focused on the relatively simpler task of analyzing the Fisher geometry of some important classes of linear models[START_REF] Hanzon | Aspects of Fisher geometry for stochastic linear systems[END_REF][START_REF] Peeters | System identification based on Riemannian geometry[END_REF]. In particular, the analytic expressions of the entries of the Fisher information matrix for several manifolds of linear MIMO systems have been discovered by Hanzon et al.[START_REF] Peeters | Symbolic computation of Fisher information matrices for parameterized state-space systems[END_REF].Let us consider first the class of stable autoregressive discrete-time processes of order 2, AR(2), in a stochastic setting in which as before the input signal is a Gaussian white noise with zero mean and unit variance. This set can be given a[START_REF] Cock | Subspace angles and distances between arma models[END_REF] In[START_REF] Belkin | Semi-supervised learning on riemannian manifolds[END_REF] in particular, dimensionality reduction is considered a factor in improving classification.

The related geometries have indeed much in common,[START_REF] Hanzon | Identifiability, recursive identification and spaces of linear dynamical systems[END_REF].

The straightforward details are neglected to improve the readability of the proof.