Xavier Pennec
email: xavier.pennec@sophia.inria.fr

Jean-Philippe Thirion

A Framework for Uncertainty and Validation of 3-D Registration Methods based on Points and Frames

In this paper, we propose and analyze several methods to estimate a rigid transformation from a set of 3-D matched points or matched frames, which are important features in geometric algorithms. We also develop tools to predict and verify the accuracy of these estimations. The theoretical contributions are: an intrinsic model of noise for transformations based on composition rather than addition; a unified formalism for the estimation of both the rigid transformation and its covariance matrix for points or frames correspondences, and a statistical validation method to verify the error estimation, which applies even when no "ground truth" is available. We analyze and demonstrate on synthetic data that our scheme is well behaved. The practical contribution of the paper is the validation of our transformation estimation method in the case of 3-D medical images, which shows that an accuracy of the registration far below the size of a voxel can be achieved, and in the case of protein substructure matching, where frame features drastically improve both selectivity and complexity.

Introduction

Many algorithms in Computer Vision concern matching tasks, whose aim is to find the correspondence between two representations of an object. Matching tasks are closely related, but not equivalent, to registration processes, which involve the evaluation of a geometric transformation. The main problem of matching methods is generally to reduce the complexity of associating features, such as in [START_REF] Ayache | HYPER : A new Approach for the Recognition and Positioning of Two-Dimensional Objects[END_REF] or [START_REF] Huttenlocher | Object Recognition using Alignment[END_REF] for Alignment, [START_REF] Grimson | Object Recognition by Computer -The role of Geometric Constraints[END_REF] for Interpretation Trees, (Lamdan [START_REF] Lamdan | Geometric Hashing : A General and Efficient Model-Based Recognition Scheme[END_REF][START_REF] Wolfson | Model-Based Recognition by Geometric Hashing[END_REF] or [START_REF] Rigoutsos | Distributed Bayesian Object Recognition[END_REF] for Geometric Hashing, [START_REF] Besl | A Method for Registration of 3-D shapes[END_REF] or [START_REF] Zhang | Iterative Point Matching for Registration of Free-Form Curves and Surfaces[END_REF] for Iterative Closest Point (ICP) methods.

In the following, we do not discuss matching methods per se, but the estimation of the 3-D motion. The traditional approach is to apply a least squares method using, for example, the singular value decomposition in [START_REF] Arun | Least-Squares Fitting of Two 3-D Point Sets[END_REF] and [START_REF] Umeyama | Least-Squares Estimation of Transformation Parameters Between Two Point Patterns[END_REF], or the quaternion representation in [START_REF] Horn | Closed Form Solutions of Absolute Orientation Using Unit Quaternions[END_REF]. Some authors, such as [START_REF] Zhuang | Robust 3-D 3-D Pose Estimation[END_REF], concentrate on the robustness of that estimation.

Uncertainty handling is a central topic in several works, like (Durrant-Whyte, 1988a;Durrant-Whyte, 1988b) or [START_REF] Ayache | Artificial Vision for Mobile robots -Stereovision and Multisensor Perception[END_REF][START_REF] Zhang | 3D Dynamic Scene Analysis: a stereo based approach[END_REF]. There are fewer studies, however, dealing with the precision of the estimated motion, which is our main concern here. Early experimental works can be found in [START_REF] Fang | Some Experiments on Estimating the 3-D Motion Parameters of a Rigid Body from Two Consecutive Image Frames[END_REF], [START_REF] Snyder | The Precision of 3-D Parameters in Correspondence-Based Techniques: The Case of Uniform Translational Motion in a Rigid Environment[END_REF], or [START_REF] Haralick | Pose Estimation from Corresponding Point Data[END_REF]. A theoretical evaluation of the errors on the rotation estimation is introduced in [START_REF] Kanatani | Geometric Computation for Machine Vision[END_REF] and [START_REF] Kanatani | Analysis of 3-D Rotation Fitting[END_REF]. The present paper is the continuation of the work presented in [START_REF] Pennec | Validation of 3D Registration Methods based on Points and Frames[END_REF].

As previous studies deal principally with point-topoint correspondences, one of our contributions is a method which applies also to frame-to-frame correspondences (a frame is a local coordinate system, i.e. a point with a trihedron). This raises some interesting problems in error modeling: a common assumption is that errors in geometric features can be modeled with an additive error on the parameters of the feature representation. We demonstrate that some important invariance properties are not conserved with this type of noise and we strongly believe that the assumption of additivity should be replaced by the composition with an error motion. This idea leads us to the devel-opment of a formalism which unifies the handling of the error for the application, composition and estimation of rigid motions, in the case of point-to-point or frame-to-frame correspondences.

Section 2 investigates rigid motions and the handling of uncertainty on geometric features, from the probabilistic point of view. Section 3 describes the statistical point of view, in particular the noise model which should be used. We present a method to estimate a rigid motion and its covariance matrix from a set of matched points or a set of matched frames. In a similar framework, we are able to estimate the fusion of several rigid motions. Those estimates are used to evaluate the noise process on the features used, and sort matched primitives (points or frames) in order to reject outliers, that is, matched primitives which are "obviously" incompatible with the estimated motion.

Section 4 describes our robust motion estimation algorithm based on the theoretical results of section 2. It can be placed after, or better, embedded within your favorite matching algorithm. Essentially, this algorithm is the iteration of three processes: motion estimation, noise process estimation and outliers rejection. It gives, from an initial set of matched primitives with or without a priori knowledge of their covariance, a final rigid motion and its covariance matrix, an estimation of the noise process on features, and the sorted set of matched primitives compatible with that motion.

Another contribution is presented in section 5, where we describe an a posteriori method to validate statistically our motion estimation algorithm. This statistical validation applies even when no "ground truth" is available.

Finally, we present in section 6 experimental results showing the performance of our registration algorithm and its statistical validation. Those results are given for synthetic data, and also for real data, namely 3D medical images, for which the evaluation of the precision of the registration is vital. Another application example is presented with 3D protein substructure matching.

Probabilistic features and rigid motions

Representation of rigid motion and frames

Let B = {o, i, j, k} be the canonical right-handed orthonormal coordinate system of Euclidean space R 3 (we have thus o = (0, 0, 0) and [i, j, k] = I d ). Let F = {t, i , j , k } be another right-handed orthonormal coordinate system (i.e. a frame), which coordi-nates are expressed in B. The rigid motion from B to F can be written in homogeneous coordinates

M = R t 0 1 with R = [i , j , k ]
Let R be a rotation of angle θ around the unit axis n. Using the rotation vector r = θn (see [START_REF] Ayache | Artificial Vision for Mobile robots -Stereovision and Multisensor Perception[END_REF] and appendix A.1), we can also represent the rigid motion M by the six dimensional vector f = (r , t ).

We note that f represents both the rigid motion M and the frame F, with respect to the canonical coordinate system B. Other representations could have been chosen for the rotation, e.g. quaternions, Euler's angles. . . but f would still represent the same rigid motion or feature frame. For convenient notation within the text, we write f = (r, t) and call it indifferently a frame or a rigid motion when there is no ambiguity.

Utility of frames in geometric computations

The use of frame features instead of points is motivated by three main reasons. Firstly, frames are natural features emerging from some geometric problems. For instance, extremal points (defined in [START_REF] Thirion | Computing the Differential Characteristics of Isointensity Surfaces[END_REF]) are points on an iso-surface optimizing a differential geometry criterion. Since they lie on a surface, they are provided with the two principal curvature directions t 1 , t 2 and the surface normal n (see Fig. 1), which forms a frame.

3-D Registration Methods

205

In protein substructure matching, an amino acid location is characterized by the position of 3 atoms in a specific geometry (see figure 1). Since the constraints are identical for every amino acid, the shape of the triangle is not discriminant and the only geometric descriptor of an amino acid is its pose, which is again described by a frame. For computer vision applications, we could also think to model each vertex of a polyhedral scene by a frame.

The second reason is that more information is available, thus higher selectivity and/or accuracy (for complexity reasons, however, the extra information should be statistically relevant in order to be useful). The introduction of normals in (Feldmar et al., 1994) allows, for instance, an easy search for seed matches by the use of bitangents, and an improvement of the accuracy in the pose estimation. In the same way, using frame features instead of points allows one to exploit the maximum information about the features (see section 6), and can lead to a drastic reduction of the problem complexity as in [START_REF] Pennec | An O(n 2 ) Algorithm for 3D Substructure Matching of Proteins[END_REF]. Last but not least, frames and rigid motions are equivalent objects from a mathematical viewpoint, and knowing how to handle uncertainty on rigid motions is also necessary.

Operations on frames and motions

Rigid motions can be composed, inverted and applied to points: the composition (operator •), inverse ( (-1) ) and apply ( ) operations on rotation vectors (see appendix A.1), can be easily extended to motions in the following way.

• Application of f = (r, t) to point x:

y = f x = r x + t • Composition of f 1 = (r 1 , t 1 ) by f 2 = (r 2 , t 2 ): f = f 2 • f 1 = r 2 • r 1 r 2 t 1 + t 2 • Inverse of f = (r, t): f (-1) = r (-1) r (-1) (-t)
Due to the equivalence of frames and motions, the application of a rigid motion f to a frame f m is simply the composition f s = f • f m . We note that in this formula, f m and f s have the same interpretation (frame or motion), but f is always a motion.

Probabilistic features

Estimated or measured values of geometric features are corrupted by measurement errors. Let x be the exact representation of such a feature and x the measured one. From a probabilistic point of view, x is the observation or realization of a n-D random vector characterized by its probability density function (pdf) ρ x . For computational reasons, however, a common assumption is to retain only the first and second order centered moments of the pdf, i.e. the expectation vector and covariance matrix of the representation considered as a random vector. This approach was mainly developed in robotics (Smith & Cheeseman, 1987;Smith et al., 1987;Durrant-Whyte, 1988a;Durrant-Whyte, 1988b) with some applications in computer vision [START_REF] Zhang | 3D Dynamic Scene Analysis: a stereo based approach[END_REF]. We note that no assumption is made about the physical noise process. The expectation and the covariance matrix of an "uncertain feature" are defined by

x = E(x) = y.ρ x (y).dy

Σ xx = E (x -x)(x -x)
The probabilistic feature x is then treated as the couple (x, Σ xx ). A common hypothesis is that the noise is centered, i.e. E(x) = x = x. In this case, a probabilistic feature is the couple (x, Σ xx ) and a deterministic feature has a null covariance matrix.

It is now interesting to see how to generalize our group operations on probabilistic features. For that purpose, we use the classical first order approximation of Jacobians.

Propagation through an explicit function Let x = (x, Σ xx ) be an m-D random vector and h a pdimensional C 1 function acting on x. Then the p-D random vector z = h(x) is determined, up to the first order, by z = (z, Σ zz ) where

z = h(x) Σ zz = J h .Σ xx .J h (1)
with

J h = ∂h(x) ∂x =    ∂h 1 ∂x 1 . . .

  

The Jacobian matrix J h generally depends on the point x where it is computed, and we denote by J h | x its value whenever the application point needs to be specified. For a function h(x, y) of two independent parameters, the mean is simply z = h(x, ȳ) and the covariance matrix

Σ zz = J h x Σ xx J h x + J h y Σ yy J h y with J h x = ∂h(x, y) ∂x and J h y = ∂h(x, y) ∂y
The composition of functions is often used to simplify the calculus. Let h 1 , h 2 and h be three multidimensional functions such as h

(x) = (h 1 • h 2 )(x) = h 1 (h 2 (x)). The Jacobian of h is then J h x = J h 1 h 2 (x) . J h 2 x
Propagation through an implicit function Let x = (x, Σ xx ) be a m-D random vector and Φ : R m × R p → R p a C 1 function. We want to investigate the p-dimensional random vector y implicitly defined by Φ(x, y) = 0. From the implicit function theorem, y = ϕ(x) exists around a given point (x 0 , y 0 ) such as Φ(x 0 , y 0 ) = 0 if and only if ∂Φ ∂y can be inverted at this point and we have in this case

∂y ∂x = ∂ϕ ∂x = - ∂Φ ∂y -1 ∂Φ ∂x
Using a first order Taylor series expansion, ȳ is implicitly defined by Φ(x, ȳ) = 0 and the covariance matrix is given by

Σ yy = ∂Φ ∂y -1 ∂Φ ∂x Σ xx ∂Φ ∂x ∂Φ ∂y - ( 2
)
Minimizing a criterion Let C be a function of class C 2 from R m × R p to R + (i.e. a criterion). We now define the p-dimensional random vector ŷ as the argument for which the criterion is minimum for a given x: ŷ = ArgMin y (C(x, y))

A necessary condition to obtain a minimum is

Φ(x, ŷ) = ∂C ∂y (x,ŷ) = 0 and H = ∂ 2 C ∂y 2 (x,ŷ)
positive definite

We are then back to the case of an implicit function and since the Hessian matrix H is symmetric, the covariance matrix of ŷ at the minimum is

Σ ŷŷ = H -1 ∂Φ ∂x Σ xx ∂Φ ∂x H -1 (3)
This propagation scheme will be used in section 2 for the rigid motion minimizing the least squares criteria.

Propagation of uncertainty for rigid motions

To use the above framework on frames and rigid motions, we need to determine the Jacobians of the three basic operations: composition, inversion of frames, and application of a motion to a point. Appendix A.1 describes the tedious computations related to rotations, in particular the Jacobian of the composition of rotation vectors.

Application of a frame (f , Σ ff ) to a point (x, Σ xx ) Let J be the Jacobian of y = f x with respect to f = (r, t) and R the rotation matrix associated to the rotation vector r.

J = ∂(f x) ∂f = ∂(r x+t) ∂(r,t) = ∂(r x) ∂r ; I 3 ∂(f x) ∂x = ∂(r x+t) ∂x = R
The computation of the Jacobian ∂(r x) ∂r is detailed in appendix A.1.4. The covariance matrix of y is then

Σ yy = J Σ ff J + R Σ xx R (4) Composition of frame (f 1 , Σ 11 ) by (f 2 , Σ 22 ) Let J 1 (resp. J 2 ) be the Jacobian of f = f 2 • f 1 with respect to f 1 = (r 1 , t 1 ) (resp. f 2 = (r 2 , t 2 )). Then J 1 = ∂(f 2 • f 1 ) ∂f 1 = ∂(r2•r1) ∂r 1 0 0 R 2 J 2 = ∂(f 2 • f 1 ) ∂f 2 = ∂(r 2 •r 1 ) ∂r2 0 ∂(r 2 t 1 ) ∂r2 I 3
This draft paper is provided to ensure timely dissemination of scholarly and technical work. Copyright and all rights therein are retained by authors or by other copyright holders. All person copying this information are expected to adhere to the terms and constraints invoked by each author's copyright. This work may not be reposted without the explicit permission of the copyright holder.

3-D Registration Methods

207

The computation of the Jacobian of r 1 • r 2 with respect to the rotation vectors r 1 and r 2 are detailed in appendix A.1.4. The covariance on f is simply given by

Σ ff = J 1 Σ 11 J 1 + J 2 Σ 22 J 2 (5)
Inversion of frame (f , Σ ff ) Let J I be the Jacobian of f (-1) = (r (-1) , r (-1) (-t)) with respect to f = (r, t). As r (-1) = -r, we have:

J I = - I 3 0 ∂(r (-1) t) ∂r (-1)
R using once again the Jacobian of the application of a rotation vector to a point (appendix A.1.4). The covariance of f (-1) is then

Σ f (-1) f (-1) = J I Σ ff J I (6)
We have now completed the computations needs to propagate the covariance matrices through the three basic operations on frames and points.

The noise model

The problem is now slightly different. In the previous sections, we considered that the pdf (or the mean value and the covariance matrix) of the probabilistic geometric features were known. In the real world, this can be a very strong assumption, in particular when we begin with simple measurements without any idea of their noise. In such a case, a natural assumption is to say that all measurements are corrupted by the same noise process. This is the notion of Identically Independently Distributed (IID) measurements [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF]. For points, the common assumption is to impose the same covariance matrix on all point measurements, which corresponds effectively to the same additive noise. Regarding frames, we have to be more careful.

Some paradoxes in standard geometric probabilities (see for instance the Bertrand paradox in [START_REF] Kendall | Geometrical probability[END_REF]) can be avoided by defining the uniform pdf over a differential manifold (i.e. the set of geometric features of a given type) as the pdf which remains invariant by the action of any transformation of a given (Lie) group. The transformations considered can be, for instance, the rigid motions: if we have a uniform distribution over our manifold, this distribution has to remain uniform when we apply the same rigid motion to all the features (for other examples, see [START_REF] Pennec | Randomness and Geometric Features in Computer Vision[END_REF]).

For our purpose, we have to impose a similar invariance property on the noise process. Let f1 and f2 be measurements of two frames f 1 and f 2 , corrupted by the same noise process. Then

f 1 = f • f1 and f 2 = f • f2 are measurements of f 1 = f • f 1 and f 2 = f • f 2 and
have to be corrupted by the same noise process (possibly different from the first one). This simply means that two identical (or comparable) distributions remain identical by a change of the reference frame, which seems reasonable. A related assumption is usually required for criterions: if y minimizes C(x, Y), it is desirable to have a criterion invariant by a given group of transformation (i.e. C(f x, f Y) = C(x, Y) for any transformation f ) in order to obtain the transformed result f y minimizing the criterion C(f x, Y) on transformed data.

Why additive noise is not adapted Let f1 and f2 be measurements of the frames f 1 and f 2 , corrupted (independently) by the same noise process of covariance Σ. If the two noises are centered, this means that fi = f i + δf i with E(δf i ) = 0 and

E δf i . δf i = Σ
If we change our view-point, i.e. apply a global motion f to exact and measured values, we obtain

f i = f • fi = f i + δf i with E(δf i ) = 0.
The first order approximation of Jacobians gives us the propagation of covariance matrices:

E δf i . δf i = J i .Σ.J i with J i = ∂(f • f i ) ∂f i
In general J i depends on f and f i , therefore J 1 = J 2 : the covariance matrix Σ 11 of f 1 is thus different from the covariance matrix Σ 22 of f 2 . The two measurements are no longer corrupted by the same noise process. As we can see, using an additive noise hypothesis can lead to paradoxes since the result of the computations depends upon the chosen reference frame.

A simple example with rotation vectors We consider an example with rotations only, represented by their ro-tation vectors. Let r 1 = ( π 2 , 0, 0) and r 2 = (0, 0, 0) be two rotation vectors with the same covariance:

Σ 11 = Σ 22 = Σ = σ 2   0 0 0 0 1 0 0 0 1  
If we apply a global clockwise rotation of π/2 around the x axis, represented by the rotation vector r = (-π 2 , 0, 0) , we obtain the new rotation vectors r 1 = (0, 0, 0) and r 2 = (-π 2 , 0, 0) . Computing the Jacobians and propagating covariances, we find that Σ 11 = 8 π 2 Σ, which is significantly different from Σ 22 = π 2 8 Σ (the ratio is approximately 1.5: see Fig. 2).

It is normal that covariance matrices change according to the motion, but since they change differently at different positions, we cannot define identical distributions as additive noises (with different means) sharing the same covariance matrix: if these matrices are identical in one coordinate system, they are not in another one. Since the IID hypothesis is used in most usual statistical methods (and will be used in the following), it is important to define it well.

A compositive noise model In order to satisfy the invariance constraint, we propose to model the measurement process by the composition of the exact value with a noise (motion): let f be the exact representation of a frame and f the measured one, then f = f •e where e is a small rigid motion around the identity. In this process, f is a deterministic frame, whereas e and f are probabilistic.

With this noise model, the invariance property is automatically verified: if g is the global motion to be applied,

f = g • f = (g • f ) • e = f • e.
We note that this compositive noise model corresponds to the left composition f • e by the exact frame. We could not have used the right composition e • f since in this case, the transformation of f = e • f by g gives f = g •e•f = e•f . We can then characterize a noise process on frames by a "random frame" e measuring the identity (i.e. the canonical frame). A centered noise corresponds to E(e) = I d , and is entirely defined by its covariance matrix: e = ( I d , Σ). From now on, the noise on frames is a compositive noise model, and we use the probabilistic measurement e of identity to represent it.

We note that this is also an additive noise in the local frame, or more specially in the tangent space of the manifold at the expected point.

Validity of the additive noise It is interesting to note that, for 2D rotations (with angular representation) or for translations (i.e. points), the composition corresponds to the addition. In these cases, compositive and additive noise models are identical, which have provided a wrong intuition for more complex cases. More generally, the additive noise model corresponds to the compositive one as soon as the action of the transformation group on features (or on itself) is linear in the considered representation, in which case the tangent space corresponds to the space itself. There exists in this case a matrix F for each motion f such that f x = F.x (or f • g = F.g) and the Jacobian ∂(f x)/∂x = F is thus independent of the position x of the feature.

The relevance of the compositive noise model to real cases will be shown in section 7 with the analysis of the estimated noise on extremal points.

Uncertainty and frames: the statistical point of view

In the first part, we focus on how to model a priori information about uncertainty, and how to propagate it through the usual geometric computations. However,

Σ 1 = Σ Σ 2 = Σ Σ 1 Σ 2
Fig. 2. The same covariance matrix on two rotation vectors in one reference frame generally leads to different covariance matrices in another reference frame. The circle is the intersection of the definition domain of the rotation vector with the (ry, rz) plane.

in real cases, we often need to estimate this information a posteriori from data. This is the goal of statistics.

The first section is devoted to the estimation of a rigid motion between two sets of matched features, while considering that the noise process on features is known. The following question is thus: how can we estimate this noise process on features? The last section of this part deals with outliers detection.

Estimation of a rigid motion

In this section, we address two main problems of estimation. The first is the estimation of a rigid motion (f , Σ ff ) from two sets of matched geometric features, namely points or frames. This is the registration problem. The second is the fusion of several measurements (f i , Σ fifi ) of the same motion in order to get a more accurate estimation. This is a data fusion problem.

Many papers on pose estimation focus on closedformed solutions in the case of point-to-point correspondences (see the references in the introduction and especially [START_REF] Kanatani | Geometric Computation for Machine Vision[END_REF]), but fewer deal with the uncertainty estimation of the pose found (see however [START_REF] Csurka | Characterizing the Uncertainty of the Fundamental Matrix[END_REF]). Moreover, we also want to estimate the rigid motion from frame-to-frame correspondences, and no classical method easily applies to such a problem.

However, a unified framework for estimation with uncertainty has been developed using the Extended Kalman Filter (EKF) in [START_REF] Ayache | Artificial Vision for Mobile robots -Stereovision and Multisensor Perception[END_REF]. We recall in section 2 the basic principles of this algorithm and introduce the notations and framework for motion estimation. In the following sections, we apply it to rigid motion estimation from point matches, from frames correspondences, and from a set of measurements of this motion. A special section (2) will be devoted to the computation of the uncertainty of the motion estimated by a standard least-squares on points.

Kalman Filtering Assume we have a set of measurements (or data) {χ i } and we search a state variable g such that, for each exact data χ i , we have the vectorial relation

z i (χ i , g) = 0
This is called the measurement equation. In our case, the state g is the sought rigid motion and the data are couples of matched points or frames, or simple measurements of this rigid motion (section 2). Since we are working with noisy data, we only know the mea-sured values of data χi = χ i + ω i (the observation). The additive noises ω i are assumed to be independent, white and centered with a known covariance Ω i :

E(ω i ) = 0 E(ω i . ω i ) = Ω i and E(ω i . ω j ) = 0 for i = j
The measurement equations are generally not linear, but assuming we know a good estimate ĝ of the state g, we can linearize them around the estimates and solve the problem with standard linear optimization techniques, namely here Kalman Filtering. This is the basis of the Extended Kalman Filtering technique. Since Kalman Filtering is a recursive filter, we assume that we have at each step i an estimate ĝi-1 of the state vector. We can then linearize the measurement equation around ( χi , ĝi-1 ) with a first order Taylor series expansion. Taking the following notations:

ẑi = z i ( χi , ĝi-1 ) M i = ∂z i ∂g = ∂z i ∂g ( χi ,ĝ i-1 ) ∂z i ∂χ = ∂z i ∂χ ( χi ,ĝ i-1 )
the Taylor expansion of the measurement equation

z i (χ i , g) = 0 gives ẑi + ∂z i ∂g .(g -ĝi-1 ) + ∂z i ∂χ .(χ i -χi ) 0
This equation can be re-written in the linear form

M i . g = γ i + ν i where γ i is the linearized measurement γ i = M i .ĝ i-1 -ẑi and ν i is a centered noise ν i = ∂z i ∂χ .( χi -χ i ) of covariance Σ ii : Σ ii = ∂z i ∂χ Ω i ∂z i ∂χ
Assuming that an initial estimate (ĝ 0 , Ĝ0 ) of the state is given, a reasonable (weighted) least-squares criterion to minimize can be constructed on the basis of Mahalanobis distances:

C = 1 2 (ĝ 0 -g) Ĝ-1 0 (ĝ 0 -g) + 1 2 i (γ i -M i .g) Σ (-1) ii (γ i -M i .g)
The recursive solution for the minimization of this criteria is called the Kalman Filter ( [START_REF] Jazwinsky | Stochastic Processes and Filtering Theory[END_REF][START_REF] Ayache | Artificial Vision for Mobile robots -Stereovision and Multisensor Perception[END_REF]). At each step, the input is an estimate (ĝ i-1 , Ĝi-1 ) of the state and the linearized measurement (γ i , Σ ii ), and the output is the updated estimation of the state (ĝ i , Ĝi ). We just recall here the equations of the filter:

K i = Ĝi-1 .M i .(Σ ii + M i . Ĝi-1 .M i ) -1 ĝi = ĝi-1 -K i .ẑ i Ĝi = ( I d -K i .M i ). Ĝi-1 (7)
For linear transformation of Gaussian variables, the Kalman Filter produces the optimal estimate of the state (in the sense of minimum error variance), which turns out to be also the maximum likelihood estimate. Moreover, it preserves the Gaussian nature of the random variables and thus there is no loss of information in keeping only the mean value and covariance matrix.

If the Gaussian assumption is removed, the Filter remains the best linear estimator, but is no longer the best one amongst all non-linear estimators. In the general case of non-linear transformations with any type of noise (which is our case in this article), the Extended Kalman Filter only represent a sub-optimal non linear estimator, but appears to provide accurate estimates in practice. However, some care has to be taken about the initial state and the order of the measurements.

Estimation from matched points Assume that we have two sets of matched points x i and y i , one transformed into the other with a rigid motion f such that y i = f x i . During the measurement process, these points are corrupted by additive noise: we only measure xi = x i + δx i and ŷi = y i + δy i and we want to estimate the motion f and the covariance Σ ff of this estimation. We assume that the errors δx i and δy i are independent with covariances Σ x i x i and Σ y i y i .

The error vector is in this case the classical difference in position z i = ŷi -f xi (and the measurement equation is z i = 0). The sought state is a rigid motion g = (r, t), and the observation vector is ( χi

, Ω i ): χi = xi ŷi Ω i = Σ x i x i 0 0 Σ y i y i
If R is the rotation matrix corresponding to the rotation vector r of the motion g, we have ∂z i /∂χ i = [-R; I d ] and hence

Σ zizi = Σ yiyi + RΣ x i x i R (8)
The EKF process is then the following.

• Initialize the state with Identity or a least squares estimate with a large covariance matrix (100 to 1000 times the least-squares covariance, for example), in order to minimize the influence of the initial state on both the estimation and its uncertainty.

• For each couple of matched points (x i , ŷi ):

Compute ĥi = ŷi -ĝi-1 xi Compute M i = ∂h ∂g = -∂(g x) ∂g estimated at (x i , ĝi-1 )
Compute Σ zizi using Eq.( 8) with rotation R i-1 . Update (ĝ i-1 , Ĝi-1 ) to (ĝ i , Ĝi ) using Eq.( 7).

Least squares In the isotropic case (more specially when Σ x i x i = Σ y i y i = σ 2 I d ), the covariance on the error vector becomes Σ z i z i = 2σ 2 I d , and the criterion reduces to a simple least squares

C = 1 2 i ŷi -f xi 2 = 1 2 i z i .z i
which can be solved exactly by several techniques. One consists in computing the barycenters of the two sets and searching for the rotation using quaternions (see [START_REF] Horn | Closed Form Solutions of Absolute Orientation Using Unit Quaternions[END_REF]).

The covariance matrix can be determined using equation ( 3) with some important simplifications: assuming that χ = (x 1 , . . . x n , y 1 , . . . y n ) is the vector of all observations and f the state, the associated implicit function characterizing an optimum is

Φ(χ, f ) = ∂C ∂f = i ∂z i ∂f z i = 0
The values of the errors z i being often small around the minimum, we can neglect in the derivatives of Φ the terms of the form z.z with respect to the terms of the form z .z (see [START_REF] Gill | Practical optimization[END_REF], section 4.7)). We obtain thus

H = ∂Φ ∂f i ∂z i ∂f ∂z i ∂f ∂Φ ∂χ i ∂z i ∂f ∂z i ∂χ
The covariance matrix of the state is then given by equation 3:

Σ ff = H -1 i ∂z i ∂f ∂z i ∂ Σ ∂z i ∂ ∂z i ∂f H -1 = H -1 i ∂z i ∂f Σz i z i ∂z i ∂f H -1
In our case, the error vectors z i are independent and have identical diagonal covariances:

z i = (0, 2σ 2 I d )
(it is under this assumption that the least-squares estimator is optimal). The central term in the previous equation collapse then into 2σ 2 H and the covariance of the motion can then be simplified in

Σ ff = 2σ 2 H -1 HH -1 = 2σ 2 H -1 (9)
with

H = i ∂( f x i ) ∂ f ∂( f x i ) ∂ f
Estimation from matched frames We now face the following problem. Assume we have two sets of matched frames {f mi } and {f si } (m for model and s for scene), transformed one into another with a global rigid motion f such that

f s i = f • f m i .
We only have access to their measured values fsi = f si • e si and fm i = f m i • e m i and we want to estimate the rigid motion f and the uncertainty Σ ff of that estimate. We assume moreover that the compositive noises are independent and centered:

e m i = ( I d , Σ m i m i ) and e s i = ( I d , Σ s i s i ).
Combining the two measurement equations with the hypothesis and isolating the error term, we get our error vector:

z i = f (-1) si • f • fmi = e (-1) si • e mi (10)
Here, we are using a particular aspect of our representation for frames: identity corresponds to a null translation and a null rotation vector. This allows us to solve this estimation problem using the EKF.

An initial estimate is easily obtained with the first couple of frames: ĝ0 = fs 0 • f (-1) m 0 . From the observation vector:

χ = fm f (-1) s Ω = Σ fmfm 0 0 Σ f (-1) s f (-1) s
and the estimations Ĵ1 , Ĵ2 , Ĵ3 and Ĵ4 at ( fmi , fsi , ĝi-1 ) of the following Jacobians:

J 1 = ∂(g•fm) ∂g J 2 = ∂(f (-1) s •(g•f m )) ∂f (-1) s J 3 = ∂(f (-1) s •(g•f m )) ∂(g•f m ) J 4 = ∂((f (-1) s •g)•f m )) ∂f m
we can simplify the Jacobians of the error vector:

∂z ∂g = J 3 .J 1 ∂z ∂χ = [J 4 ; J 2 ]
and its covariance matrix is given by

Σ zz = Ĵ4 Σ f m f m Ĵ 4 + Ĵ2 Σ f (-1) s f (-1) s Ĵ 2 (11)
The EKF estimation can thus be summarized as:

• Initialize the state with

(ĝ 0 , Ĝ0 ) = ( fs 0 • f (-1) m0 , Σ (fs 0 •f (-1) m 0 )(fs 0 •f (-1) m 0 ) )
• For each couple of matched frames ( fm i , fs i ):

Compute ẑi = f (-1) si • ĝi-1 • fm i and the Ja- cobians Ĵ1 , Ĵ2 , Ĵ3 , Ĵ4 . Compute M i = Ĵ3 . Ĵ1 and Σ zizi using Eq.(11). Update (ĝ i-1 , Ĝi-1 ) to (ĝ i , Ĝi ) using Eqs.(7).
Fusion of rigid motions or frames Consider now a set { fi } of measurements of the same frame f such as fi = f • e i . The goal is to estimate the frame (f , Σ ff ). Assuming an error e i = ( I d , Σ ii ) on each measurement, the error vector is

z i = f (-1) i • f = e (-1) i (12)
An initial estimate is easily obtained with the first frame: ĝ0 = f0 . The measurement χ is now simply χ = f (-1) and its covariance matrix

Ω = Σ f (-1) f (-1) .
The EKF process can thus be summarized as:

• Initialize the state with (ĝ 0 , Ĝ0 ) = ( f0 , Σ f 0 f 0 ) • For each measurement of frame fi : Compute ẑi = f (-1) i • ĝi-1 Compute M i = ∂(f (-1) •g) ∂g and Ĵ1 = ∂(f (-1) •g) ∂f (-1) estimated at ( f (-1) i , ĝi-1 ) Compute Σ z i z i = Ĵ1 Σ f (-1) f (-1) Ĵ 1 Update (ĝ i-1 , Ĝi-1 ) to (ĝ i , Ĝi ) using Eqs.(7).

Estimation of the noise process on features

In the previous section, we assumed that the noise process corrupting the measurements was known, and we estimated the motion. We are now trying to estimate the noise process while the motion is known (or already estimated). In fact, we only estimate its covariance matrix since the mean value should be null with our hypotheses (which is verified experimentally).

Estimation of the noise on points

Anisotropic noise: Let {x i } and {y i } be two sets of matched points with measured values xi = x i + δx i and ŷi = y i + δy i , and f be an estimate of the motion f = (r, t) linking them together (y i = f x i ).

The error vector z i = ŷi -f xi is, up to the first order, z i = δy i -r δx i if we assume an exact motion (this assumption will be removed later).

Assuming that each measurement error δx i and δy i comes from noise processes of covariances Σ xx and Σ yy , the covariance on the error vector z i should be Σ zz = Σ yy + RΣ xx R , which can be estimated by

Σzz = 1 N i z i z i (13)
If we assume that one set of points is exact, then we can solve for the uncertainty of the other:

Σ xx = 0 =⇒ Σyy = Σzz and Σ yy = 0 =⇒ Σxx = R Σzz R
If we assume that both sets have the same noise Σ, we obtain Σzz = Σ+RΣR , which is uniquely solvable for Σ unless R is a rotation of angle θ = π/2 or θ = π (see [START_REF] Koch | Parameter Estimation and Hypothesis Testing in Linear Models[END_REF] for a solving method). If the rotation between the two sets is small, the covariance on points is approximated at first order by

Σxx = Σyy = Σzz 2 = 1 2 N i z i z i
Isotropic noise: We have in this case Σ xx = Σ yy = Σ = σ 2 I d , and hence the covariance matrix of the error vector should be Σ zz = 2 σ 2 I d . Using the estimation proposed in equation ( 13) and taking the trace, we obtain:

Tr( Σzz ) = 2σ 2 . d = 1 N z i 2
where d is the dimension of the space (d = 3 in our case). Since Ĉ = 1 2 z i 2 is the value of the leastsquares criterion at the minimum, the variance estimation can be summarized as

σ2 = Ĉ d N

Estimation of the noise on frame features

Standard noise process: With the same notations as in the previous sections, and neglecting once again the error on the estimate f of the motion f , the error is:

e i = f (-1) si • f • fm i = e (-1) si • e m i
Assuming a common underlying process of measurement errors e s i = e m i = ( I d , Σ), we can determine that e (-1) s i = ( I d , Σ) (since the Jacobian of the inversion is -I d at the origin), hence e i = ( I d , 2Σ). It should be noted that the covariances behave here as usual, but only because errors are close to identity. An estimator of the covariance Σ = Σ m i m i = Σ s i s i of frame features is given by

Σ = 1 2 N i e i e i
A simplified model: In some cases, for example with a small number of matches, the estimation of the above covariance matrix can be unstable. We can use instead a kind of isotropic model of noise with a standard deviation σ θ for the rotation part and σ d for the translation part. Splitting the error vector into a rotation and a translation component e i = (e θi , e di ), we can then estimate the variances by

σ2 θ = e θ i 2 6 N and σ2 d = e d i 2 6 N
This draft paper is provided to ensure timely dissemination of scholarly and technical work. Copyright and all rights therein are retained by authors or by other copyright holders. All person copying this information are expected to adhere to the terms and constraints invoked by each author's copyright. This work may not be reposted without the explicit permission of the copyright holder.

Discussion about noise estimation: Several problems arise in the covariance estimation scheme: assuming a Gaussian noise model with outliers (contaminated Gaussian for instance), it is clear that estimating the covariance matrix on all points, including outliers, will produce a significant elevation of the covariance.

On the other side, estimating on inliers only will induce an under-estimation which corresponds, on a theoretical point of view, to the following partial integration on the inlier domain D

Σ = D (x -x).(x -x) .ρ(x).dx
instead of the whole space R d for the real covariance matrix Σ.

With the Gaussian assumption, and a χ 2 threshold on the Mahalanobis distance selection (see the following section), this equation reduces to Σ = α(χ).Σ where the correction factor α is defined by an integral of the χ 2 distribution. In practice, this factor is only slightly inferior to one with standard values of the χ 2 threshold, and we prefer to compute the covariance with matches that have a Mahalanobis distance less than 1.5 or 2 times the χ 2 threshold. This includes only the outliers that are almost inliers and ensures a very small underestimating bias with usual values of the χ 2 threshold. This remains, however, the weak point of our method concerning robustness, at least theoretically (see section 4.2).

Another bias in the covariance estimation comes from the lest-squares technique we choose to determine the motion (which is indeed not exact): if one set of features is an exact model, the residuals are a good estimation of the real error, but when both sets of features are corrupted by noise, the residue minimization makes the residuals slightly inferior to the real error. From [START_REF] Bard | Nonlinear Parameter Estimation[END_REF], the number of observations N in all our covariance estimations should thus be replaced by N = N -l/m, where m is the dimension of the vectorial equations and l the number of parameters we have estimated. In our case, this is a rigid motion estimation, so l = 6. If we estimate it from point matches, we have m = 3 whereas m = 6 from frames. We should then replace N with N is the previous equations, where

• N = N -2 for points, • N = N -1 for frames.

Rejecting outliers: Mahalanobis distance and χ 2 test

Another interesting problem that arises within this framework is compatibility. For a given motion, we want to know, for instance, if a scene feature is compatible with a model feature, i.e. if the scene feature can be considered as the transformation of the model feature modulo measurement errors. Considering sets of matched features, we may also want to sort the matches by saliency. These two problems are important for the estimation since least squares techniques are known to be sensitive to outliers, and the EKF is moreover order dependent: a preliminary phase where matches are sorted by relevancy and outliers rejected often gives substantial improvements on the quality of the estimation.

A related question is to determine when two measurements originate from the same object and differ only because of measurement error, This can be used, for instance, for clustering motions, i.e. to find subgroups of consistent motions. A fusion could then be used within each subgroup to obtain a better estimate of the motion.

Such problems are traditionally tackled via Mahalanobis distances for sorting hypotheses and χ 2 tests for rejecting outliers. Until now, we did not need to assume Gaussian distributions on our random variables. For the χ 2 test, however, this assumption is required. Two main reasons allow us to assume that the distributions of the noise is Gaussian: firstly, the measurement errors are often the sum of independent errors and thus tend toward a Gaussian process by the Central Limit Theorem. Secondly, when the mean and variance of an unknown distribution are the only information available (which is our case), a simple maximum entropy derivation gives the Gaussian distribution as the one that assumes the minimal information [START_REF] Bard | Nonlinear Parameter Estimation[END_REF].

Two points and a motion Let (x, Σ xx ) and (y, Σ yy ) be two measured points and (f , Σ ff ) a rigid motion. We want to test the validity of the hypothesis y = f x.

We saw in section 2.4 how to compute (f x, Σ (f x)(f x) ) from (x, Σ xx ) and (f , Σ ff ). If the motion was independent of the data x and y, then the covariance on the error vector z = y -f x should be zero with covariance Σ = Σ yy + Σ (f x)(f x) . In our case, the motion is computed from the data and the independence assumption does not hold. Thus we have to use the covariance of the residue [START_REF] Förstner | Reliability Analysis of Parameter Estimation in Linear Models with Applications to Mensuration Problems in Computer Vision[END_REF]:

Σ = Σ yy + ∂(f x) ∂x .Σ xx ∂(f x) ∂x
which does not take into account the uncertainty due to the motion (it is often very small with respect to the features uncertainty). A χ 2 test is then well suited to verify if z is compatible with its theoretical covariance.

The squared Mahalanobis distance µ 2 is

µ 2 = z Σ (-1) z
The statistical test says that the hypothesis is good if µ 2 ≤ ε, where the threshold ε is set such that, if the hypothesis is true, we will choose it with the probability α. Some values of ε and α for this 3-D χ 2 test are given found in table 1: for example, if µ 2 > 11.34, there is less than 1% of probability that y be a measurement of f x. We can thus reject this match.

Two frames and a motion Let (f m , Σ mm ) and (f s , Σ ss ) be two measured frames and (f , Σ) a rigid motion. We want to test the hypothesis that f s = f •f m . From section 2.4 we can compute the covariance Σ on the error vector e = f (-1) s

• f • f m (with the same remark as for points concerning the independence of the motion) and use the squared Mahalanobis distance µ 2 = e Σ (-1) e and a 6-D χ 2 test to decide if the hypothesis is true or not (see table 1).

Two frames Let (f 1 , Σ 11 ) and (f 2 , Σ 22 ) be two measured frames. We want to test the hypothesis that they are measurements of the same frame f . Let Σ be the covariance on the error vector e = f (-1) 1

• f 2 . We can once again use the squared Mahalanobis distance µ 2 = e Σ (-1) e and a 6-D χ 2 .

Conclusion: uncertainty and rigid motions

We introduce in this section a new model of noise based on composition and show why it is better adapted than the classical additive noise. Estimation problems within this framework require the computation of a covariance matrix of the estimate: the Extended Kalman Filter provides a unified formalism to handle this. At This framework can be used in a large number of vision problems. We present in the next section an application to registration; but we could also use this framework for matching algorithms (see [START_REF] Pennec | An O(n 2 ) Algorithm for 3D Substructure Matching of Proteins[END_REF]). Further applications include the clustering of rigid motions between image structures, or rigid motions in a displacement field, in order to detect and isolate objects. This can also prove to be useful when dealing with the accumulation stage of geometric hashing or Hough transform algorithms.

A practical and robust algorithm for motion estimation

The output of a feature-based matching algorithm usually consists in two sets of matched features. The basic idea is that if we can compute a reliable estimation of the motion, and the confidence that we have in it, we can give a confidence value on the registration at each point of the image. The method is the following. We take as input of the evaluation step two sets of matched features and compute the rigid motion between them. Assuming that all the features are corrupted by the same noise process, we compute statistics (i.e. the covariance matrix of the process) on features, and use this information to recompute the rigid motion and its associated covariance. This iterative process can be continued until convergence. This allows us to predict a variance on the position of each registered point of the image.

The iterative process

The considered features are again points and frames, but initially without any uncertainty information 1 . Without loss of generality, such a framework can be extended to other types of features. In order to use an EKF, we first need to estimate the features covariance matrix (this is detailed in section 3.2). We start with 1. A least squares estimate of f , using point positions only, or directly a least median of squares on frames (see section 4.2). 2. Estimate the noise process on features. 3. Order matches by increasing Mahalanobis distances.

We now have enough information on the features to start the iterative process. Step 1 can be ignored if a reliable motion is provided, as step 3 if matches are already sorted with respect to some other criterion. This ordering of matches is, however, needed to ensure the robustness of the EKF: since the equations are linearized in the EKF, the estimate is order-dependent. If the differences are usually not significant, some special cases can lead to quantitative errors. Step 2 is not used if an estimation of the noise process (even a gross estimation) is given on input (which is what we usually do unless it is the first time we use a new type of data).

Ordering also helps outliers rejection for the motion computation while keeping these matches for other purposes. In protein substructure matching, for instance in [START_REF] Pennec | An O(n 2 ) Algorithm for 3D Substructure Matching of Proteins[END_REF], we want to compute the motion with reliable amino-acids matches (which can have some precise interactions with other molecules), but we also want to see less constrained matches that point out other structural similarities.

The iterative process itself is as follows.

1. Estimate the motion (f , Σ ff ) from the matches using an EKF (in decreasing significance order) until the Mahalanobis distance of the matches becomes larger than a predetermined ε threshold (this is the χ 2 test to reject outliers). 2. Estimate the noise process with all features (see section 3.2). 3. Sort matches by their Mahalanobis distance. Outliers are the tail of the list.

The process is repeated until convergence, or for a preset maximum number of iterations (typically 5 to 10). Convergence means in particular that the current estimate of the motion is exactly the same as at the previous step.

The estimation of the noise process is the weak point of our scheme for robustness since it is relatively errorprone. This step is dismissed if a reliable estimation of this noise is furnished, but when it comes to a new type of data, we need to estimate it. In practice, we estimate the noise on features (i.e. its covariance matrix) only once per three or four iterations in order to stabilize the convergence (in this case around 15 to 20 iterations), and results on real data appear to be quite accurate. A rigorous treatment of this point would require the robust estimation of a covariance matrix.

Variations on the algorithm

Since a frame is composed of a point and a trihedron, we can also apply point techniques to frames. On this basis, we can then distinguish four ways of estimating both the motion and the noise process: EKF on frames with complete or simplified noise model, EKF on points with a complete noise model (on points), and standard least-squares on points (isotropic noise model).

The appropriate method depends upon the number of matches, the relative quality of the trihedron and the isotropy of the noise. A comparison of the different methods is presented in section 6.1.

In terms of robustness (see [START_REF] Huber | Robust Statistics[END_REF][START_REF] Meer | Robust Regression Methods for Computer Vision: A Review[END_REF] for a review on robust statistics), the algorithm we propose belongs to the class of iterative redescending M-estimators (i.e. robust algorithms based on least squares). However, the breakdown point (the maximum amount of outliers before failure) is not very high. The algorithm can be robustified to a breakdown point of about 0.5 by using a least median of squares between the initialization step (in order to have an estimation of the noise process) and the iterative process, that achieve the goal of efficiency (i.e. reaches a variance on the estimate close to the lowest possible one). This is particularly easy and cheap in complexity for frames. Indeed, a single match determines a unique motion between the model and scene sets: for each of the N frame matches, we classify the other matches by increasing Mahalanobis distances and take the value of the median for the score of this match. The match minimizing the median squared Mahalanobis distance is then used as the initial motion for the iterative process. The complexity of O(N 2 log N ) can be drastically reduced by a Monte Carlo sampling (see [START_REF] Meer | Robust Regression Methods for Computer Vision: A Review[END_REF]).

Final Precision evaluation

From ( f , Σ ff ), we can compute for each point x of the model image (or object) its transformation ŷ = f x and the uncertainty Σ ŷŷ = Σ ( f x)( f x) . This information could be used directly as an input for other statistical algorithms but is too rich to give a simple idea of the level of accuracy for the end user. In order to give a more intuitive information, we compute instead the RMS error expectation at this point: if y = f x is the real location of the transformed point x, we shall find it in ŷ = f x. The expectation of the squared distance d 2 = y -ŷ 2 between them is

σ 2 (ŷ) = E (ŷ -y) (ŷ -y) = Tr(Σ ŷŷ )
In the case of 3-D medical images, and in order to characterize the precision with a single value, we can estimate (among other measurements) the averaged RMS for the 8 vertices of the 3D image (we call this value the typical boundary precision σ corner ), which is representative of the maximum error, or estimate the averaged RMS on the position of the matches frames (the typical object precision σ object ).

Conclusion on registration and accuracy

From two sets of matched features, and a χ 2 threshold ε, we propose a method to estimate the motion f and its precision Σ ff , discard outliers, and provide an estimation of the noise process on features. Then we extract from the motion uncertainty a single value characterizing the registration accuracy.

A statistical validation of registration methods

We have now a way to estimate motions and associated covariances, but what confidence can we have in it ? We propose in this section a method to estimate the accuracy of a registration algorithm when no "ground truth" is available. It also validates the estimation process for feature based registration methods which generate covariance matrices. This is a statistical method which only assumes that the noise on features is centered.

External markers

We believe that no method can give an exact reference for the motion, but when a method A is one order of magnitude more accurate than a method B, A can serve as a reference for an experimental validation of B. This means that all methods are ultimately statistical.

Debates about the use of external markers are typical in the case of 3D medical image registration problems. Here, validation methods are vital (literally), because the registration can be used to plan a surgical procedure. In fact, external markers have a very strong advantage in matching methods, because we can extract features corresponding to physical and discriminable objects and obtain an "exact" matching. But regarding registration, we have to keep in mind that the measurements of marker positions are corrupted just as the measurement of any other landmark, which implies the elaboration of a theoretical model of the markers, and the application of a statistical method to evaluate the motion and its uncertainty. Other authors have concluded experimentally that, for the case of high resolution medical images, marker-based techniques are probably much less accurate than image-based techniques (see [START_REF] Van Den Elsen | Multimodality Matching of Brain Images[END_REF]).

Hence, marker-based techniques are not more deterministic than other registration methods, and the results of the present paper could also be applied to them (generally, markers are modeled as points or frames).

Analytical methods

Ideally, we might think of modeling analytically the whole process, from the physical object positions to the final estimated motion. For real cases, this implies modeling object deformations (nothing is really rigid), evaluating the distortions due to image acquisition and reconstruction process, evaluating the errors in the extraction of feature points, before considering the errors made during registration. This is generally impossible to do in practice, when all this also depends on the shape of the object, and the tuning of the acquisition device.

An a posteriori estimation of the errors

The method that we present now applies when there is no ground truth for the motion, and when the analytical determination of the errors is impossible (the majority of cases). We consider the registration method as a black box which takes two representations of the same object as input, and gives an estimated rigid motion as output.

Assuming independence of couples of matched features, we can subsample the matches in order to get several independent estimates of the unknown motion, and compare them. In particular, we can split the set of matches in two sets approximately equal in size and obtain two estimates of f :

f 1 = f • e 1 and f 2 = f • e 2
The two motions should be very close, and we can study their "difference" e = ( f 2 ) (-1) • f 1 = (e 2 ) (-1) • e 1 which does not depend on the exact motion f and should be close to identity.

If matches are randomly subsampled to a fixed number M of matches, in order to conserve independence and a similar distribution of matches in space, then the estimation errors e 1 and e 2 follow the same law of covariance Σ(M ). If moreover both the feature extraction and the registration method are unbiased 2 (which means that E(e 1 ) = E(e 2 ) = I d ), then we can derive that the covariance of (e 2 ) (-1) is also Σ and the measured error e should then have a covariance Σ ee (M ) = 2Σ(M ).

Repeating the experiment with n homogeneous sets of matches, we can estimate the covariance matrix of the estimation of f from M feature matches by

Σ(M ) = 1 2 n n i=1
e i e i

A posteriori validation of the uncertainty

Our estimation scheme produces not only a motion estimation f but also a covariance estimation on it Σ ff . The goal of this section is to verify a posteriori that the distribution of f has for mean value f with the given covariance. We consider first the case of synthetic data where the exact motion f is known. Since every experiment is based on a different exact motion and produce a different covariance matrix on the estimate, we have to "normalize" our results in order to obtain several estimates of the same distribution.

To normalize the mean motion, we compute for each registration the error vector e = f (-1) • f (for simpler notations, the index of the registration experiment is omitted). Its distribution should have a zero mean, but its covariance matrix is given, since f is exact, by Σ ee = J.Σ ff .J with J = ∂e/∂ f . A new change of variables is needed to normalize the distributions with respect to this covariance: under the Gaussian hypothesis, the Mahalanobis distance with identity µ 2 = e Σ (-1) ee e should be χ 2 6 distributed. We can now repeat this experiment on M pairs of images to obtain M independent values µ 2 i and verify if it is really χ 2 6 distributed. The Kolmogorov-Smirnov test [START_REF] Press | Numerical Recipices in C[END_REF] is well adapted to do that (referred from now on to the K-S test), but since it only gives a binary answer, we also use the fact that the mean value of a χ 2 6 distribution is 6 and its variance is 12: We call validation index the estimated mean value of µ 2 i :

I = μ2 = 1 M µ 2 i
and its variance is computed accordingly with

σ 2 I = 1 M -1 (µ 2 i -μ2 ) 2
This index can be interpreted as an indications on how the estimation method under-estimates (I > 6) or over-estimates (I < 6) the error on the estimated motion. It is a kind of relative error on the error estimation.

This method can be generalized to validate our estimation scheme with real data: randomly splitting the set of matches in two (approximately equal) sets, we compute two estimates ( f 1 , Σ f 1 f 1 ) and ( f 2 , Σ f 2 f 2 ) which are two independent measures of the exact motion f . Their "difference" e = ( f 2 ) (-1) • f 1 can now be used just as above. Some results with synthetic and real data are given in sections 6.1 and 6.2.

Conclusion on statistical validation

We show in this section that all validation methods are ultimately statistical and present a quite general method to estimate the quality of a registration process, independently of our previous method. This leads us to the validation of our uncertainty prediction with both synthetic and real data.

Experiments

The first part of this section deals with the comparison of the four different methods on synthetic data, and shows how to choose the best method among isotropic or anisotropic noise and points or frames features. We will see through these experiments that all the methods are perfectly validated in their own domain. The second and third part focuses on real data, with 3D medical image registration and protein substructure matching.

Synthetic data experiments

For our synthetic experiments, we have used lists of associated frames or points with a Gaussian error distribution on the positions and on the orientations and a uniform random placement in 256x256x256 voxels This draft paper is provided to ensure timely dissemination of scholarly and technical work. Copyright and all rights therein are retained by authors or by other copyright holders. All person copying this information are expected to adhere to the terms and constraints invoked by each author's copyright. This work may not be reposted without the explicit permission of the copyright holder. images (one voxel equal one millimeter). The exact motion is chosen randomly but the translation is limited to a third of the image size in order to obtain an overlap between the two images. The covariance matrix of the simulated noise process is given by the experiments on 3D medical images (see the covariance matrix of table 2). We use a χ 2 value of 16 for frames and 8 for points, which corresponds to a 96% confidence interval, and limit the number of iterations to a maximum of 15.

We also run several other sets of experiments using a real model randomly moved and perturbed, and a uniform model of noise with the same covariance. All those experiments showed very similar results.

Anisotropic versus simplified noise model

We have presented in section 3.2 two models of noise for points (iso and anisotropic) and their equivalent for frames. The question we investigate here is which one should we use, and does the anisotropic model provide an improved accuracy of the registration? In these experiments, frames are generated with a compositive model of noise and points with an additive one.

We claim that the relevant parameter to decide which method is the most suited is the number of matches. Indeed, the estimation of the 21 parameters (resp. 6 for points) of the covariance matrix is less stable than the estimation of the 2 (resp. 1) parameters of the simplified model, in particular with a small number of samples. However, for a sufficient number of matches, we can estimate the whole covariance matrix with a good confidence, even if the original noise is isotropic.

In figures 3 and 4, we plot on the top the validation index (ideally 6) and its standard deviation. The dotted lines represent the theoretical mean value and standard deviation of a χ 2 6 distribution. We also indicate experiments where the K-S test at 5% rejects the validation. On the bottom of the figures, we plot the typical precisions as a function of the number of matches for isotropic methods and the gain in precision obtained using anisotropic methods. Each sample point is the averaged result of 100 trials.

As expected, the validation index shows that isotropic methods need less matches to have a good uncertainty prediction than anisotropic methods ones: about 10 points are needed for the isotropic method on points to be reliable, versus 15 to 20 matches for the anisotropic one. For frames, these figures are as expected a little bigger: respectively 15 to 20 and 30 to 40 matches are needed. On the other hand, anisotropic methods lead to a more accurate motion estimation: in this case, it is roughly 1.25 times more accurate for frames and between 1.15 and 1.4 times more accurate for points. However, we note that these values are just examples and do vary with the anisotropy of the original noise. From now on, we use the simplified model of noise for less than 40 matches with frames (20 with points) and the anisotropic one otherwise.

The difficult and possibly inaccurate point in our our scheme is therefore the covariance matrix estimation which is, as we have already noted, not really robust. In order to verify the robustness and accuracy of the remainder of the algorithm, we run another series of experiments where we fix the covariance on features (i.e. we consider that it is already known and we do not have to estimate it) and introduce about 10% of outliers. In all cases, even with 5 points or frames matches, the Kolmogorov-Smirnov test (at 5%) successfully validate the rigid motion and its uncertainty estimation.

Comparison of Frame and Point based methods An interesting question is what improvement in accuracy is

brought by the use of frames instead of points ? Since we add trihedra to points in order to make frames, frame methods should be at least as accurate as point ones. On the other hand, the extraction of trihedra from real images is usually based on higher derivatives than the point position. Hence, points are usually less noisy than trihedra and we do not expect a breakthrough for accuracy in using frames. However, as points become noisier and noisier, there should be sensible improvement with frames.

We choose to run the experiment with a random number of points between 150 and 250, a fixed standard deviation of σ θ = 0.02 rad for the orientation sampling, and a standard deviation increasing from σ d = 0.2 mm to 2.8 mm on point position sampling.

Since the typical precisions do vary with the number of matches and the noise process, but in a similar way for both methods, we plot in figure 5 the ratio of typical precisions (standard deviations) for point and frame with respect to the standard deviation on position. Each sample is the averaged result of a set of 50 synthetic "images".

Using frames in spite of points can then lead to an sensible improvement in accuracy, particularly with a small number of matches, or when points becomes noisy with respect to trihedrons. On the other hand, several experiments showed that the frame based method remain at least as accurate as the point based one even when trihedrons become very noise. In general, we should then use the maximal amount of information we have on features and thus the frame based method if we can define trihedrons.

Table 2. Estimated covariance matrix for the compositive noise on extremal points (expressed in the local frame).

Σ ee = e r e t          
0.0024 0.0000 -0.0000 0.0002 -0.0011 0.0000 0.0000 0.0030 -0.0000 0.0001 -0.0000 0.0000 -0.0000 -0.0000 0.0373 -0.0006 -0.0003 0.0001 0.0002 0.0001 -0.0006 0.2276 0.0011 0.0008 -0.0011 -0.0000 -0.0003 0.0011 0.3157 -0.0015 0.0000 0.0000 0.0001 0.0008 -0.0015 0.0838

          er e t
This draft paper is provided to ensure timely dissemination of scholarly and technical work. Copyright and all rights therein are retained by authors or by other copyright holders. All person copying this information are expected to adhere to the terms and constraints invoked by each author's copyright. This work may not be reposted without the explicit permission of the copyright holder.

Real data experiments: 3D medical imaging

We present results from an experiment performed using 3D Magnetic Resonance images (MRI). The images were provided by R. Kikinis from the Brigham and Woman's Hospital, and are part of an extensive study of the evolution of the Multiple Sclerosis (MS) disease. The same patient recieved a complete 3D MR examination several times during one year (typically 24 different 3D acquisitions). The aim is to register Fig. 6. The same slices from four registered 3D MR images of the Multiple Sclerosis study. Note the evolution of two MS lesions (white spots). One is growing in the anterior left hemisphere, one is shrinking in the right posterior hemisphere. There is two weeks between each acquisition.

precisely in 3D all those images in order to segment the lesions and evaluate very accurately their evolution (MS lesions are white spots in the images in figure 6).

The images are first echo, 256 × 256 × 54 voxels, the voxel size is 1 mm × 1 mm × 3 mm. We have already presented an algorithm to perform, fully automatically, the registration of the images in [START_REF] Thirion | Extremal Points : definition and application to 3D image registration[END_REF]. However, there is a need to estimate accurately the errors in these registrations. The registered images are resampled using a tri-linear interpolation method. The same slices of several different 3D MRI are presented in figure 6, after resampling with the estimated motions. Figure 7 presents the differences between images and shows visually the quality of the registration, and also that the evolutions of the lesions are clearly detected.

Points and frames in medical images Our registration algorithm relies on the extraction of feature points in 3D medical images, defined with differential geometry criteria (see figure 8). In our case, these are the Extremal Points, as defined in [START_REF] Thirion | The 3D Marching Lines Algorithm : new results and proofs[END_REF], which are those points of the object surface for which both principal curvatures are extremal. The interesting thing is that not only do we get some invariant measurements associated with those points (the principal curvatures), which are used to reduce the complexity of the matching process, but we get also the principal Fig. 7. Differences between images after registration, with respect to the first one. The intensity is multiplied by 5, and shifted such that no difference is a gray value. The growing lesion appear as a white disk, and the shrinking lesion as black a disk. directions, which form, with the normal to the surface and the extremal point itself, an orthonormal basis, that is, a frame (see figure 1 in section 2.2).

Typically, we extract about 2000 extremal points from a 3.5 million voxels image. Our matching algorithm produces about 600 pairs of associated extremal points between two images, with a residual mean square error (RMS) of about 1mm. With a density fewer than 0.1% of extremal points in the image, the probability of false matches is very low.

Results

With 24 different images, we could have generated 576 couples of images. In order to achieve independence, we have restricted our experiments to 58 randomly chosen registrations. For each registration, we randomly split the list of matches in two approximately equal lists and compute the validation index I from the two rigid motion estimates. We fuse then the two motion estimates in order to compute the final motion and its uncertainty.

We find a typical boundary precision around σ corners = 0.114 mm and a typical object precision far below the voxel size: σ object = 0.060. The validation index is I = μ = 6.19 with a variance of σ 2 I = 17.05 (remember that the theoretical values are 6 and 12) and the K-S test validate our results for real data with a significance level of 0.16. These values were found using frames and showed an improvement of 10% in accuracy with respect to points.

Analysis of the estimated noise on features These experiments showed an interesting difference between the point and frame noise model. Indeed, point based methods gave a quite isotropic noise of standard deviation σ = 0.5, whereas the frame based method gave the covariance matrix on frame features displayed in Table 2 .

As said in the synthetic experiments section, this covariance is roughly diagonal, with standard deviations σ r x = 0.05, σ r y = 0.055, σ r z = 0.20 for the rotation vector part (in radian), and σ d x = 0.5, σ d y = 0.55, σ dz = 0.25 for the position part (in mm), which gives a mean standard deviation on position of σ = 0.46, comparable with the one computed on points only.

In order to interpret these values, we have to remember that, with the compositive noise model, the error vector is expressed in the local frame. In our case, the frame is defined by {t 1 , t 2 , n} where t 1 and t 2 are the two principal directions of the surface at the extremal point and n the normal. As far as the rotation vector is concerned, the values of the standard deviations show that the expected error rotation around the n (= z) axis is four times the expected error rotation around the t 1 or t 2 axes. This means that the normal n is about four times more stable than the principal directions. Regarding the position, we can see that the coordinate along the normal to the surface is 3 to 4 times more stable than coordinates in the tangent plane, which is in accordance with what was expected for extremal points.

The compositive noise model exhibits here an anisotropy which is not detectable with a classical noise: the absence of orientation information leads to a quasi isotropic noise for points, and an additive noise on frames would fail to detect this regularity since the error would not be expressed in local frames. This effect therefore constitutes an a posteriori justification of our compositive noise model with real data.

Discussion Our matching algorithm tries to estimate the motion of a single rigid substructure, with the largest number of common extremal points, which is the brain in those images. However, with that level of accuracy, the assumption of a global rigid motion does not hold any more, and we can distinguish several structures undergoing slightly different motions. The skull, for example, can move (quite rigidly) with respect to the brain, and the skin surface is subject to large deformations (this is very clear with animated sequences spanning the year study). The points corresponding to those structures have been rejected as outliers, but, with the help of the framework developed in this paper, we can think to classify and discriminate substructures. This is an improvement of our registration method that we are currently implementing.

Even for the brain itself, there are local deformations, for example due to the studied MS lesions, which are larger than the typical boundary error we measured, but because we are interested in an average motion for the whole structure of interest (here the brain), it makes sense to search for such an accurate motion. In fact, our registration is so accurate that it allows us to visualize, for the first time, the dynamic effect of the lesions on the surrounding brain tissues.

3D substructure matching of proteins

Most biological actions of proteins depend on some typical parts of their three-dimensional structure, called 3D motifs. To automatically discover corresponding 3D motifs between proteins, we proposed in [START_REF] Pennec | An O(n 2 ) Algorithm for 3D Substructure Matching of Proteins[END_REF] a new 3D substructure matching algorithm based on the geometric hashing technique. The key feature of the method is the introduction of a 3D reference frame attached to each amino acid (see figure 1 in section 2.2). This allows to compute for every couple of amino acids 6 invariants, that are used in the first step of the algorithm: a geometric hashing aimed to recognize possible correspondences between individual model and scene frames. The introduction of frame features instead of points drastically reduce the complexity of both the preprocessing and recognition stages, typically from O(n 4 ) to O(n 2 ). A clustering step based on motions performs then the aggregation of compatible matches. Individual motions within a cluster are fused using the scheme of section 2 in order to obtain a unique and accurate global motion. The last step is an alignment test realized as a kind of "iterative closest frame" (by analogy with ICP [START_REF] Besl | A Method for Registration of 3-D shapes[END_REF]), where our robust motion estimation scheme is embedded.

The noise on frame features is fixed on input and determines the type of common substructures the algorithm extracts. A small noise produces small but very accurate sets of matches, whereas a large noise favors more global substructure matches. For all our experiments, we use the atom coordinates of proteins provided by Brookhaven National Laboratory's Protein Data Bank [START_REF] Bernstein | The Protein Data Bank: A Computerbased Archival File for Macromolecular Structures[END_REF][START_REF] Abola | Protein Data Bank[END_REF]. Visualization is done using the RasMol program of R. Sayle [START_REF] Sayle | RasMol: A program for fast realistic rendering of molecular structures with shadows[END_REF].

Detection of the Helix-Turn-Helix motif Structural motifs can be defined as the super-secondary structure. They are the simple combination of a few secondary structure elements. Some of them are associated with particular functions or are simple parts of larger structural and functional assemblies. For instance, the Helix-Turn-Helix motif is responsible for the binding of DNA within many procaryotic proteins. Some of them bind tightly to the DNA at a promoter of a gene, preventing RNA polymerase from fixing and hence blocking the initiation of the transcription. They are repressors. Conversely, activators bind next to the promoter and help polymerase to bind.

We choose to compare two proteins known to bind DNA: the tryptophan repressor of E. Coli (PDB code 2WRP [START_REF] Lawson | Flexibility of the DNA-Binding Domains of TRP Repressor[END_REF], 105 amino-acids) and phage 434 CRO (PDB code 2CRO [START_REF] Mondragon | Structure of Phage 434 Cro Protein at 2.35 Angstroms Resolution[END_REF], 65 amino-acids), whose Helix-Turn-Helix sequence were determined in [START_REF] Brennan | The Helix-Turn-Helix DNA binding motif[END_REF][START_REF] Harrison | DNA recognition by proteins with the Helix-Turn-Helix motif[END_REF].

Looking for a 3D binding motif, we used a quite small model of isotropic noise (σ θ = 0.1 rad = 5 deg and σ d = 0.35 Å ). The algorithm ends up with only the correct matches from (15 MET -66 MET) to (36 ALA -87 ASN). The last match (37 GLY -88 SER) is indeed very arguable considering the distance after registration and especially the difference in orientation. The typical object precision due to the registration (on the 22 matched amino-acids) is given to σ obj = 0.29 Å. We show in figure 9 the two proteins and the registration found in figure 10. The only two other common substructures found score 13 and 8 matches and correspond to alpha helices, which are very stable secondary structure elements.

Discussion In order to test the stability of our algorithm, we also did the experiment with a noise two times larger (σ θ = 10 deg and σ d = 0.7 Å ). We just find four additional matches preceding the beginning of the HTH motif ( from (7 LYS -61 LEU) to (11 ILE -64 GLY)). The two other clusters now score 14 and 9 matches. This shows that the detection of the HTH motif is very stable, and we argue that this is mainly due to the selectivity brought by the use of frames.

Indeed, the orientation of an amino acid is crucial to determine the position of collateral chains and most interactions of the protein happen within these side atoms. The position of these atoms is then not only determined by the position of the backbone but also its orientation and using just points to represent amino acids generally leads to a significant amount of additional matches with non compatible orientations. This implies a drastic reduction of selectivity for the matching process. In this case, frames bring rather more selectivity than more precision.

Conclusion

We developed a new formalism to handle uncertainty for rigid motion evaluation, which consider residual errors as residual geometric motions, and not as an additive noise applied to the measurements coordinates or Fig. 10. Registration found between the HTH motifs of 2CRO and 2WRP. We can see that not only the backbone is very well matched, but also collateral chains are pretty well conserved.

to the parameters of the unknowns. This leads to new formulae for the equations of measurement which are used within the Extended Kalman Filter framework to evaluate both the motion and its covariance matrix, in the case of 3-D points or 3-D frames matches. We provide a quite general scheme to estimate approximately the accuracy of a registration method, and a more precise way to validate the predictions of our registration scheme. Several experiments on both synthetic and real data validate our framework and justify the compositive noise model. A practical result is to show that, in the case of 3D medical images, a precision of the registration far below the voxel size can be achieved. In the case of 3D protein substructure matching, frame features drastically improve the selectivity and reduce the complexity of the process.

This work also opens new theoretical questions about the noise on geometric features under a given type of motion, and offer numerous possibilities of extensions: other types of features, more general transformation groups, and also other high level algorithms based on this framework, such as clustering or matching.

two representations are given by Rodrigues formula (see for instance 1993)).

R = I d + sin θ S n + (1 -cos θ)S n 2 = cos θ I d + sin θ S n + (1 -cos θ)n.n (A1)
The matrix S n is the skew matrix corresponding to (left) cross product: for all vector v we have S n v = n×v. If the coordinates of n are (n x , n y , n z ), the matrix S n is

S n =   0 -n z n y n z 0 -n x -n y n x 0   and we have the relation S n 2 = n.n -I d , used
to derive the second part of equation A1. Note that S n uniquely determines the vector n. Conversely, let Tr(R) be the trace of R; the parameters are:

     θ = arccos Tr(R) -1 2 S n = R -R 2 sin θ (A2)
The last equation is valid only when θ ∈]0; π[. Indeed, for θ = 0 (i.e. identity) the rotation axis n is not determined, and sin(θ) = 0 for reflections (θ = π).

R close to identity: θ is small Since the axis n is not defined for identity, there is a singularity and a numerical instability around it. However, we can compute the rotation vector with a Taylor expansion:

S r = θ S n = θ 2 sin θ (R -R ) = 1 2 1 + θ 2 6 (R -R ) + O θ 4
R close to a reflection: π -θ is small The axis is this time well defined, but we have to use another equation.

From Rodrigues formula, we get

R + R -2 I d = 2(1 -cos θ)S n 2 , and since S n 2 = n.n -I d , we have n.n = I d + 1 2(1 -cos θ) R + R -2 I d Let = (1 -cos θ) -1
; taking diagonal terms in the last equation gives n 2 i = 1 + .(R i,i -1) and thus

n i = ε i 1 + .(R i,i -1)
The off diagonal terms are used to determine the signs

ε i : considering that the sign of n 1 is ε ∈ {-1; +1}, we can compute that sign(n k ) = ε.sign(R 1,k + R k,1 )
If we have an exact reflection, the sign ε does not matter since rotating clockwise or counter-clockwise gives the same result, but for a quasi-reflection, this sign is important. In this case, the vector w = 2 sin θ n is very small but not identically null: it can be computed without numerical instabilities with S w = R -R . Since θ < π, the largest component w k in absolute value of this vector must have the same sign as the corresponding component n k of vector n.

A.1.2. Rotation vector

The representation of a rotation by its axis and angle values or by a unit quaternion are very useful for a lot of problems, but they are not minimal: there are quadratic constraints needed to enforce unitary vectors. This is a problem for most minimization techniques, and especially the extended Kalman filter (EKF).

A minimal representation is obtained with the rotation vector r = θ n. Some problems of uniqueness are encountered around θ = π and four charts are theoretically needed to define an atlas (a set of differentiable maps covering the complete manifold). The rotation vector can be computed from the rotation matrix with the equations of section A.1.1. Conversely, Rodrigues formula allows us to compute the rotation matrix from the rotation vector:

R = I d + sin θ θ S r + (1 -cos θ) θ 2 S r 2
with θ = r . In order to avoid numerical instabilities, we have to use a Taylor expansion for a small θ:

sin θ θ = 1 -1 6 θ 2 + O θ 4 (1-cos θ) θ 2 = 1 2 -1 24 θ 2 + O θ 4
A more developed presentation of the rotation vector is given in [START_REF] Ayache | Artificial Vision for Mobile robots -Stereovision and Multisensor Perception[END_REF].

Operations on rotation vectors Let r be a generic rotation vector. The associated rotation matrix is denoted by R and the by θ = r and n = r/ r ). We are interested in this section in applying r to the vector x: y = r x, inverse it: r (-1) , and compose rotation vectors r 1 and r 2 : r = r 2 • r 1 . The two first operations are easy to perform:

r x = R.x and r (-1) = -r
The composition is more complex. Of course we can compute the associated rotation matrices, multiply them, (R = R 2 .R 1 ) and come back from R to the rotation vector r, but it would be difficult to differentiate. We choose to use the unit quaternion representation as an intermediate step (see section (A.1.3) below). Let r 1 and r 2 be two rotation vectors. We compute the associated unit quaternions q 1 and q 2 , and multiply them: q = q 2 * q 1 . Then we come back to the rotation vector r(q): r 2 • r 1 = r(q(r 2 ) * q(r 1 ))

A.1.3. Quaternions

Quaternions are the elements of a 4 dimensional algebra on R which we shall note Q. This is also the first non-commutative division ring (skew field) found by Hamilton in 1843. We can construct this algebra in several ways, but the we are interested in considers a quaternion q ∈ Q as a couple q = (a, v), where a ∈ R is the real part and v ∈ R 3 the so-called pure part. The operations defined on quaternions to form the algebra are :

• Addition : (a 1 , v 1 )+(a 2 , v 2 ) = (a 1 +a 2 , v 1 +v 2 )

• Internal multiplication :

(a 1 , v 1 ) * (a 2 , v 2 ) = (a 1 a 2 -< v 1 |v 2 >, v 1 × v 2 + a 1 v 2 + a 2 v 1 )
where '×' and '< .|. >' are the usual cross and dot products on R 3

Moreover, we define the conjugate quaternion and the norm

• (a, v) = (a, -v) • |q| 2 = q 2 Q = q * q = a 2 + v 2 R 3 = q 2 R 4
This allows to write very simply the inverse quaternion q -1 = q/|q| 2 for q = 0. We shall note that the norm is compatible with the product:

|q 1 * q 2 | = |q 1 |.|q 2 |.
The set of quaternions (0, x) with x ∈ R 3 is trivially identified with R 3 itself. Let x and y be two vectors (elements of R 3 ). Their quaternion product is x * y = (-< x|y >, x × y) and q * x * q is a vector for any quaternion q.

A more detailed introduction to quaternions and their properties is available in [START_REF] Casteljau | Les quaternions[END_REF] and [START_REF] Altmann | Rotations, Quaternions, and Double Groups[END_REF].

Quaternions and rotations Let q be a unit quaternion. Then there exists θ ∈ [0, π] and n unit vector on R such that q = cos θ 2 , sin θ 2 n . The map

R q : Q -→ Q x -→ y = q * x * q
is an inner automorphism of Q that conserves pure quaternions (null real part). Its restriction to R 3 is the vectorial rotation of R 3 with angle θ around the unit vector n. In a symmetric way, we can match to every rotation of R 3 two unit quaternions q and -q.

Let R ∼ q denote the association between rotation matrix R and rotation quaternion q. As direct properties of this representation, we have:

• If R 1 ∼ q 1 and R 2 ∼ q 2 then R 1 .R 2 ∼ q 1 * q 2 . • If R ∼ q then R -1 ∼ q
By definition, the application of R ∼ q to vector x is

• y = Rx = q * x * q

Conversions between unit quaternions and rotation vectors Let r = θ n and q = ±(a, v) be related to the same rotation R. We have then a = cos(θ/2) and v = sin(θ/2) n. The conversion from rotation vector to quaternion is then q(r) = cos(θ/2) ; sin(θ/2) θ r with θ = r . Conversely, the rotation vector can be written (the notation q =1 = = = means equal if q = 1) r(q) = 2 sign(a) arcsin v

a 2 + v 2 v v q =1 = = = 2 sign(a) arcsin( v ) v v
More details will be found with Jacobians of these operations.

A.1.4. Jacobians

We now have to differentiate the previous operations. Jacobians (see section 2.4) are needed to propagate the covariance matrices, and in the linearization of the Extended Kalman Filtering equations. We recall first the differentiation of some usual operators for 3D vectors.

• Norm and Normalization: Application of a rotation vector: r x = R.x The problem is to compute the Jacobian ∂(r x) ∂r .

∂ x ∂x = x x ∂ ∂x x x = -S x
Let α, β, γ and δ be the following functions of θ with their Taylor expansions for a small θ:

α = sin θ/θ = 1 -θ 2 6 + O θ 4 β = (1 -cos θ)/θ 2 = 1 2 -θ 2 24 + O θ 4 γ = α /θ = (cos θ -α)/θ 2 = 1 3 -θ 2 30 + O θ 4 δ = β /θ = (α -2β)/θ 2 = -1 12 + θ 2 180 + O θ 4
With Rodrigues formula, we can then write

Composition of two rotation vectors

The most complex part of these computations is the following: how to compute the Jacobians of r 2 •r 1 with respect to the two rotation vectors? From r = r 2 •r 1 = r(q(r 2 ) * q(r 1 )), we can easily derive with the chain rule ∂r ∂r 1 = ∂r ∂q . ∂q ∂q 1 . ∂q 1 ∂r 1

(A6)

The formula being symmetrical for the derivation with respect to r 2 . The problem is now to compute these intermediate Jacobians.

From rotation vector r to unit quaternion q. Let q = (a, v) be one of the unit quaternion associated with From unit quaternion q to rotation vector r. Let q = (a, v) be a quaternion. The unit quaternion q q represent the rotation associated with rotation vector r. We can get r from q with the following equation.

r(q) = 2sign(a) arcsin v

a 2 + v 2 v v q =1 = = = 2 sign(a) arcsin( v ) v v
If a is positive, we can derive that

∂r ∂a = -2v a 2 + v 2 q =1 = = = -2v ∂r ∂v = 2a a 2 + v 2 v.v v 2 - 2 v arcsin v a 2 + v 2 S v 2 v 2 q =1 = = = 2 arcsin( v ) v I d +2 a - arcsin( v ) v v.v v 2
For a negative, the above equations can be used for q = -q and thus Composition of the two quaternions: q = q 2 * q 1 Let q i = (a i , v i ). We have then

q = q 2 * q 1 = a 1 a 2 -< v 1 |v 2 > v 2 × v 1 + a 2 v 1 + a 1 v 2
The Jacobians can then be easily written

∂q ∂q 1 = ∂(q 2 * q 1 ) ∂q 1 = a 2 I d +   0 -v 2 v 2 S v2   (A9)
and (beware of the minus sign for the right-bottom block of the matrix due to the cross product):

∂q ∂q 2 = ∂(q 2 * q 1 ) ∂q 2 = a 1 I d +   0 -v 1 v 1 -S v 1   (A10)
This draft paper is provided to ensure timely dissemination of scholarly and technical work. Copyright and all rights therein are retained by authors or by other copyright holders. All person copying this information are expected to adhere to the terms and constraints invoked by each author's copyright. This work may not be reposted without the explicit permission of the copyright holder.

Fig. 1 .

 1 Fig. 1. Frame features naturally arise from geometric problems. Top: (extremal) points on iso surfaces. Bottom: geometric description of an amino acid.

Fig. 3 .

 3 Fig. 3. Validation index with respect to the number of matches for anisotropic (top right) and simplified (top left) noise model on point features. The typical precisions are presented for the isotropic model (bottom left), along with the ratio between isotropic and anisotropic precisions (bottom right).

Fig. 4 .

 4 Fig. 4. Validation index with respect to the number of matches for anisotropic (top right) and simplified (top left) noise model on frame features. The typical precisions are presented for the isotropic model (bottom left), along with the ratio between isotropic and anisotropic precisions (bottom right).

Fig. 5 .

 5 Fig. 5. Ratio between Point and Frame based precisions methods with respect to standard deviation on position. The error bar is the standard deviation in 50 trials.

Fig. 8 .

 8 Fig. 8. Lines of extremal curvature on the brain circumvolutions.extremal points, which are used in the registration, are specific points of those lines.

Fig. 9 .

 9 Fig. 9. The CRO (2CRO) of phage 434 on the left and the tryptophan repressor of E. Coli (2WRP) on the right. The matched part (the HTH motif) is displayed in black.

  Inversion of a rotation vector: r (-1) = -r The differentiation of the inversion is very simple since:

  r x = x + α.S r .x + β.S r 2 .xand taking into account the following derivatives∂(Sr.x) ∂r = -S x ∂(Sr 2 .x) ∂r = S x .S r -2S r .S xand ∂θ/∂r = r /θ, we can differentiate it using the chain rule, which gives after factorizations∂(r x) ∂r = -S x γ.r.r -β.S r + α. I d -S r .S x δ.r.r + 2β. I d (A4)Link with other works: Another way to differentiate the rotation application is given by N. Ayache in (Ayache, 1991): let η be the function η = (1 -α)/θ 2 (η = 1 6 -θ 2 120 + O θ 4 for a small θ), U r the matrixU r = η S r 2 + β S r + I dand du the infinitesimal vector du = U r .dr. Then he showed that, for any small increment dr of the rotation vector r, the increment dR of the rotation matrix R is given by dR = S du .R. To differentiate the application, we can then write(R + dR)x -R.x = dR.x = S du .R = -S (R.x) .du = -S (R.x) .U r .drand hence∂(r x) ∂r = -S (R.x) .U r (A5)which is only a factorized form of equation (A4). From a computational point of view, the first form is much faster, in particular if we note that S r .S x = x.r -< x|r > I d = x.r -Tr(x.r ) I d

  rotation vector r, and θ = r . From a = cos(θ/2)

  To summarize the calculi for q = 1, let µ = v and τ and υ the following functions of a and µτ = 2 sign(a) v) = -2v ; τ I d + υ.v.v (A8)

Table 1 .

 1 Table of the χ 2 distribution for 3 and 6 degrees of freedom.

	α Dim 50%	90 %	95 %	99%
	ε	3	2.37	6.25	7.81	11.34
	ε	6	5.35	10.65 12.59	22.46

This draft paper is provided to ensure timely dissemination of scholarly and technical work. Copyright and all rights therein are retained by authors or by other copyright holders. All person copying this information are expected to adhere to the terms and constraints invoked by each author's copyright. This work may not be reposted without the explicit permission of the copyright holder.

Acknowledgements

We would like to thank Nicholas Ayache, Mike Brady Zhengyou Zang for stimulating discussions about error handling, and also Dr Ron Kikinis from the B&W Hospital, Harward Medical School, Boston, who provided us with the MR images. Part of that research has been supported by the Esprit Basic Research Action VIVA, and also by a contract with Digital Equipment.

Notes

1. If the covariance matrix (the noise process) is already provided with the feature, we can use it directly. 2. The absence of bias for the registration process can be verified on synthetic data.

Appendix

A.1. Rotations

A.1.1. Geometric parameters: axis and angle

Let R be a 3D rotation matrix (RR = I d and det(R) = 1). It is characterized by its axis n (unit vector) and its angle θ. The relationship between these This draft paper is provided to ensure timely dissemination of scholarly and technical work. Copyright and all rights therein are retained by authors or by other copyright holders. All person copying this information are expected to adhere to the terms and constraints invoked by each author's copyright. This work may not be reposted without the explicit permission of the copyright holder.

Pennec and Thirion

Computations concerning rotations are now completed: we can apply, inverse and compose rotation vectors and compute the Jacobians of these operations.


Alexandre Guimond

Gérard Subsol

Jean-Philippe Thirion

Automatic MRI Database Exploration and Applications

published or not. The documents may come   L'archive ouverte pluridisciplinaire