Xavier Pennec
email: xpennec@sophia.inria.fr

Toward a Generic Framework for Recognition Based on Uncertain Geometric Features

Keywords: 3-D object recognition, invariants of 3-D objects, feature-matching, uncertain geometric features, generic features, matching error analysis

The recognition problem is probably one of the most studied in computer vision. However, most techniques were developed on point features and were not explicitly designed to cope with uncertainty in measurements.

The aim of this paper is to express recognition algorithms in terms of uncertain geometric features (such as points, lines, oriented points, or frames). In the first part we review the principal matching algorithms and adapt them to work with generic geometric features. Then we analyze some noise models on geometric features for recognition, and we consider how to cope with this uncertainty in the matching algorithms. We then identify four key problems for the implementation of these algorithms. Last but not least, we present a new statistical analysis of the probability of false positives that demonstrates a drastic improvement in confidence and complexity that we can obtain by using geometric features more complex than points.

Introduction

The recognition problem is probably one of the most studied in computer vision. (See for instance [START_REF] Besl | Three-dimensional object recognition[END_REF][START_REF] Chin | Model-based recognition in robot vision[END_REF].) Many algorithms have been developed to compare two images or to recognize objects against an a priori model. These models are mostly constructed from point features and can vary according to rigid or affine transformations. One can cite for instance the Geometric Hashing [START_REF] Lamdan | Geometric hashing : A general and efficient model-based recognition scheme[END_REF], the ICP [START_REF] Besl | A method for registration of 3-D shapes[END_REF][START_REF] Zhang | Iterative point matching for registration of free-form curves and surfaces[END_REF], or the alignment algorithm [START_REF] Ayache | Hyper : A new approach for the recognition and positioning of two-dimensional objects[END_REF][START_REF] Huttenlocher | Object recognition using alignment[END_REF].

However, while these algorithms are efficient for well-scattered data with little noise, it becomes more and more difficult to find small objects in a complex and noisy scene. The maximal complexity of these algorithms is then reached and the problem of uncertainty handling becomes critical. Moreover, the geometric models of the real world often lead one to consider features that are more complex than points, such as lines [START_REF] Grimson | On the sensitivity of geometric hashing[END_REF], planes [START_REF] Faugeras | Three-Dimensional Computer Vision: A Geometric Viewpoint[END_REF], oriented points, or frames [START_REF] Pennec | A geometric algorithm to find small but highly similar 3D substructures in proteins[END_REF][START_REF] Pennec | A framework for uncertainty and validation of 3D registration methods based on points and frames[END_REF]. As we will show below, using these kinds of features directly in the recognition algorithm can lead to important improvements in complexity, accuracy, and robustness. On the other hand, we must be very rigorous in handling uncertainty to avoid false negatives and paradoxes [START_REF] Pennec | Uniform distribution, distance and expectation problems for geometric features processing[END_REF].

The aim of this paper is express the above recognition algorithms in terms of uncertain geometric features. We restrict our analysis to the matching of similar features in the same space (3D-3D or 2D-2D for instance). In the following section, we briefly describe the theory already developed for handling the uncertainty of geometric features. In Section 3, we review the principal matching algorithms and adapt them to work with generic geometric features. In Section 4, we analyze three noise models for recognition (bounded, probabilistic, and a combination of both) and how to modify the matching algorithms to cope with uncertainty in features measurement. This leads to the identification of four key problems for the implementation in Section 5. In the last section, we analyze the drawbacks of uncertainty on recognition algorithm: false positives. We present a new method to evaluate qualitatively their probability of occurrence and show that using geometric features more complex than points drastically improves the algorithms' performance and robustness.

Geometric Features

We have shown in [START_REF] Pennec | Uniform distribution, distance and expectation problems for geometric features processing[END_REF] that geometric features generally do not belong to a vector space but rather to a manifold and that this induces paradoxes if we try to use the standard techniques for points with them. For instance, we can represent a 3-D rotation by its matrix R, the corresponding unit quaternions ±q, or the rotation vector r = θ • n. Using the barycenter to compute the mean, we obtain either R =1 n i R i , q = 1 n i q i , or r = 1 n i r i . The three results correspond to different rotations. 1 From a mathematical point of view, all the operations we define on features should rely on intrinsic characteristics of the manifold and not on the vector properties of some particular chart. Moreover, there generally is a transformation group acting on the manifold that models the possible image viewpoints and/or the subject movement in the images. Any operation should be invariant or "covariant" with respect to the action of this group. For instance, the barycenter of points is invariant under the action of rigid or affine transformations: A( x i ) = Ax i . In the case of geometric features, these requirements are much more difficult to meet, and we have to look at rather deep results in differential geometry in order to understand the intrinsic characteristics of the manifold we can use.

In this section, we summarize the essence of the theory developed in [START_REF] Pennec | L'Incertitude dans les Problèmes de Reconnaissance et de Recalage: Applications en Imagerie Médicale[END_REF]. More specifically, we present some intrinsic characteristics of the manifolds of geometric features and we show that most useful operations on such features can be expressed in this framework using a very small number of operations (atomic operations). The interested reader can consult [39, chap. 9], [START_REF] Klingenberg | Riemannian Geometry[END_REF], and [START_REF] Do Carmo | Riemannian Geometry. Mathematics[END_REF] for more theoretical results.

Riemannian Manifolds

In the geometric framework, one specifies the structure of a manifold M by a Riemannian metric. This is a continuous collection of dot products on the tangent space at each point x of the manifold. Thus, if we consider a curve on the manifold, we can compute at each point its instantaneous speed vector and its norm, the instantaneous speed. To compute the length of the curve, we can proceed as usual by integrating this value along the curve. The distance between two points of a connected Riemannian manifold is the minimum length among the curves joining these points. The curves realizing this minimum for any two points of the manifold are called geodesics. The calculus of variations shows that geodesics are the solutions of a system of second-order differential equations depending on the Riemannian metric.

In this article, we assume that the manifold is geodesically complete, i.e., that the definition domain of all geodesics can be extended to R. This means that the manifold has no boundary nor any singular points that we can reach in a finite time. As an important consequence, the Hopf-Rinow-De Rham theorem states that there always exists at least one minimizing geodesic between any two points of the manifold (i.e., whose length is the distance between the two points).

Exponential chart

From the theory of second-order differential equations, we know that there exists one and only one geodesic starting at a given feature x with a given tangent vector. This allows us to develop the manifold in the tangent space along the geodesics (think of rolling a sphere along its tangent plane at a given point). The geodesics going through this point are transformed into straight lines, and the distance along these geodesics are conserved (at least in a neighborhood of x). The tangent planes at points x and y of the sphere S 2 are different: the vectors v and w of T x M cannot be compared to the vectors t and u of T y M. Thus, it is natural to define the dot product on each tangent plane. Right: The geodesics starting at x are straight lines in the exponential map, and the distance along them is conserved. The function that maps to each vector the corresponding point on the manifold is called the exponential map. This map is defined in the whole tangent space T x M (since the manifold is geodesically complete), but it is one-to-one only locally around the origin. If we look for the maximal domain where this map is oneto-one, we find out that it is a star-shaped domain delimited by a continuous curve C x called the tangential cut-locus. The image of C x by the exponential map is the cut locus C x of point x. Roughly, this is the set of points where several minimizing geodesics starting from x meet. 2 On the sphere S 2 for instance, the cut locus of a point x is its antipodal point, and the tangential cut locus is the circle of radius π .

The exponential map within this domain realizes a chart called the exponential chart. It covers the whole manifold except for the cut locus of the development point, which has a null measure. Let -→ xy be the representation of y in this chart. Then its distance to the origin is dist(x, y) = y x . This chart is somehow the "most linear" chart of the manifold with respect to the feature x: geodesics starting from this point are straight lines, and the distance is (locally) conserved along them.

Homogeneous Manifolds

Invariant metric Now, since we are working with features onto which acts a transformation group that models the possible image viewpoints, it is natural to choose an invariant Riemannian metric. (The existence conditions are detailed in [START_REF] Pennec | L'Incertitude dans les Problèmes de Reconnaissance et de Recalage: Applications en Imagerie Médicale[END_REF].) This way, all the measurements based on distance are independent of the image reference frame or the transformation of the image: dist(x, y) = dist(g x, g y).

Let o be a point of the manifold that we call the origin: we call principal chart the exponential chart at the origin for the invariant metric and we denote by x the representation of x in this chart. Assuming that we have chosen an orthonormal coordinate system, the distance with the origin is dist(o, x) = x . Now, let f y be a "placement function," i.e., a transformation such that f y o = y. The distance between any two points is dist(x, y) = dist f (-1) y x, o = f (-1) y x . One can verify that this formula is independent of the choice of the placement function.

Exponential chart and map at other points

The previous formula shows that the vector (f (-1) y x) is the representation of x in an exponential chart at feature y. The coordinate system of this chart is orthonormal, but it depends on the chosen placement function. It is often more interesting to use the nonorthogonal coordinate system induced by 2. More precisely, the cut locus is the closure of this set.

VIDERE 1:2

Toward a Generic Framework for Recognition 60 the principal chart, which is independent of the placement function: this is the vector

-→ yx = J (f y ) • (f (-1) y x), where J (f x ) = ∂( f x) ∂ x x=o
is the Jacobian of the origin's translation in the principal chart. From a practical point of view, this means that we can reinterpret the local calculations on points as local calculations in the principal chart of our manifold by replacing xy with (-1) • S δx ). For instance, the empirical covariance matrix with respect to a point y, becomes [START_REF] Pennec | L'Incertitude dans les Problèmes de Reconnaissance et de Recalage: Applications en Imagerie Médicale[END_REF]:

-→ yx and x + -→ δx with exp x ( -→ δx) = f x (J (f x )
x x ( y) = 1 n n i=1 -→ yx i • -→ yx i T = 1 n J (f y ) •   n i=1 f (-1) y x i • f (-1) y x i T   • J (f y ) T .
(1)

Atomic operations

In fact, we have shown that most of the interesting operations for us on deterministic and uncertain geometric features can be expressed in the principal chart of the manifold with only the few operations we have introduced. From a computational point of view, this means that the only operations we have to implement for each type of feature are the action of a transformation (f x) with its Jacobians

∂(f x) ∂f
and ∂(f f) ∂ x , and the placement function f x with its Jacobian ∂(f x ) ∂ x . Every higher-level operation can be expressed as a combination of these and thus can be implemented independently of the considered type of feature. To simplify further computations, we can add to these atomic operations the Jacobian of the translation of the origin J (f x ).

Practical Implementation

Computer model of random features Let x be a random feature. To define the mean value, we have to replace the standard expectation by the Fréchet expectation [START_REF] Pennec | Uniform distribution, distance and expectation problems for geometric features processing[END_REF]:

E [ x ] = arg min y∈M E dist(y, x) 2 .
Assuming that this mean value is unique, let E x = x be its representation in the principal chart. Its covariance matrix in the same chart is given by

xx = E -→ xx • -→ xx T = J ( f x) • E ( f (-1) x x) • ( f (-1) x x) T • J ( f x) T .
This information is often sufficient to model the random feature and provide a compact representation. Thus, from a computational point of view, we define a random feature by its approximation: x ∼ ( x, xx ). In this framework, a deterministic feature has a null covariance matrix. Random transformations are modeled similarly.

Basic operations

Using this representation, the action of a random transformation f ∼ ( f, ff ) on a (random) feature x ∼ ( x, xx ) gives the random feature

y = f x ∼ f • x, yy , VIDERE 1:2
Toward a Generic Framework for Recognition 61

where yy = J f • ff • J T f + J x • xx • J T x with J f = ∂( f x) ∂ f f= f and J x = ∂( f x) ∂ x x= x.
The distance between two features is simply

dist(x, y) = -→ xy x = f (-1) x y T • f (-1)
x y .

Similarly, the Mahalanobis distance between a random feature x ∼ (x, xx ) and a deterministic feature y becomes

µ 2 (x, y) = -→ xy T • (-1) xx • -→ xy = (f (-1) x y) T • J (f x) T • (-1) xx • J (f x) • (f (-1)
x y).

The Mahalanobis distance between two random features is defined by

µ 2 (x, y) = min z µ 2 (x, z) + µ 2 (y, z) .
Higher-level operations Based on the atomic and basic operations, one can define many higher-level operations on geometric features. For instance, we developed in [START_REF] Pennec | Computing the mean of geometric features: Application to the mean rotation[END_REF] three gradient descent algorithms to compute the mean feature by minimizing the standard and Riemannian least-squares distance, the weighted Riemannian least-squares distance, and the Mahalanobis distance. We also develop in [START_REF] Pennec | Registration of uncertain geometric features: Estimating the pose and its accuracy[END_REF] similar algorithms for computing the optimal registration.

Summary of Geometric Features

In short, the important point to note is that geometric features usually belong to manifolds that are not vector spaces. However, developing the manifold along the geodesics onto its tangent space at the origin give a chart (the principal chart) which is almost linear: geodesics going through the origin are straight lines, and the distances are conserved along them. Using an invariant metric for the transformation group acting on features, we can translate this chart at any point of the manifold. From a practical point of view, this means that we can reinterpret the local calculations on points as local calculations in the principal chart of our manifold by replacing

-→ yx = x -y with -→ yx = J (f y ) • (f (-1) y
x), and

x + -→ δx with exp x ( -→ δx) = f x (J (f (-1)
x) ) • δx. It turns out that all the interesting operations on deterministic and probabilistic features and transformations can be expressed in the principal chart using only a few atomic operations and their Jacobians: the composition and inversion of transformations (f • g and f (-1) ), the action of a transformation (f x), and the placement function (f x ). On this basis, we can construct many higher-level algorithms, such as the computation of the mean feature or the computation of the transformation between two sets of matched features (registration).

Example of features

In [START_REF] Pennec | L'Incertitude dans les Problèmes de Reconnaissance et de Recalage: Applications en Imagerie Médicale[END_REF], we have implemented this framework for 3-D rigid transformations acting on different kinds of features such as frames, semi-oriented frames, and points. Frames are composed of a point and an orthonormal trihedron and are equivalent to rigid transformations. The principal chart is composed of the rotation vector representing the trihedron (or the rotation) and a vector representing the point position (or the translation). Semi-oriented frames model the differential properties of a point on a surface. They are composed of a point

Matching Algorithms

In this paper, we consider that matching is aimed at finding correspondences between two set of features, whereas registration is the computation of the transformation between two matched sets of features. These two problems are often intrinsically linked to constitute the feature-based recognition problem. Here, we focus on matching methods.

We have one set X = {x 1 , . . . x m } ∈ M m of m features modeling the first image and similarly a set Y = {y 1 , . . . y n } ∈ M n of n features modeling the second image. To distinguish the two sets, we use the vocabulary from the model-based recognition: we call model the set X, as if it was the model of an object stored in a library, and scene the set Y, in which we are looking for objects.

If we had two ideal acquisitions of the same object (with also perfect low-level processing), the model and the scene would have the same number of features and the two sets would be simply transformed into one another using a transformation f ∈ G, with possibly a permutation in index of features. The matching problem is in this case very simple: find the index permutation that identifies features up to a transformation f. In practice, this is not so easy: the object we are looking for is usually not the only one in our images, and hence there are numerous extra features in the two images that do not have to be matched but that we cannot suppress a priori. This is called clutter. Moreover, even if an object is perfectly known (from its CAD model for instance), an image can be noisy enough to inhibit the detection of some features. We call this occlusion as in computer vision, even if the phenomenon is different in other domains such as 3-D imaging. Last but not least, the measure of features is inherently noisy and the superimposition of matched features after registration will never be perfect.

In general, one considers that an object is recognized if there exists a transformation f ∈ G that superimposes (within a given error) a sufficient number of features in the two images. Thus, we want to maximize both the number of matches and the quality of these matches, that is the goodness of fit after registration.

Let π be the matching function that associates to the index i of a model feature x i the index j = π(i) of the matched feature y j in the scene or the null feature * if it is not matched. This function can thus be considered as an application from {1 . . . m} to {1 . . . n + 1} or from X to Y ∪ { * }. There are (n + 1) m such functions in the correspondence space . Some constraints can be added, such as the symmetry of matches. The basic idea in matching is to maximize the number of matches, but matches should also be consistent with the configuration of features. One often uses the following simple plausibility function to estimate the likelihood of a match: α(z, z ) = 1 if dist(z, z ) ≤ ε and 0 otherwise. The matching function is thus defined as the one that maximizes the matching score (f x is the action of the transformation f on the feature x):

π = arg max π∈   i α f x i , y π(i)   . ( 2
)
At the same time, we want to optimize the quality of matches and find the transformation that superimposes the model onto the scene. This is often done using least squares, with the convention that dist(x, * ) = 0.

f = arg min f∈G   i dist f x i , y π(i)   . (3)
There is thus a joint search in the correspondence space (find π ∈ ) and in the transformation space (find f ∈ G). Each of these problems taken independently is relatively simple, but solving both together is much harder. In the remainder of this section, we investigate in the sequel the main algorithms used in image processing to solve this problem and how to formulate them in terms of features.

Interpretation Trees

If we look at the equations, we find that searching for the transformation f is linked to the matching solution π . The interpretation trees method proposed in [START_REF] Grimson | Object Recognition by Computer: The Role of Geometric Constraints[END_REF] relies on the dual method: looking for π in the correspondence space and computing a posteriori the transformation f to validate the interpretation. The basic algorithm is thus to assume at each node of the interpretation tree the match of an unused feature in X. When we arrive at a leaf we validate or discard this interpretation by looking for the transformation.

Several classical methods in artificial intelligence allow one to prune the search in this exponential tree ((n + 1) m leaf). One of the most interesting is the introduction of geometric constraints with unary invariants. Assume that our features are segments: the length of a segment in the scene can only be inferior (modulo the error) to the length of the model segment since the only modifications on the length are due to occlusion and measurement errors: this is a unary constraint on the matches. Now if two segment matches satisfy the unary constraints, the angle between these segments in the model and in the scene should be almost identical: this is a binary constraint of invariance. For each type of feature, we can similarly develop a set of unary, binary, and possibly higher-order geometric constraints that will prune the interpretation tree and guarantee a local consistency during the depth-first search. The search of the transformation will actually give the proof of the global consistency on leaves. The introduction of measurement error is quite simple in this scheme: we just have to propagate the error on features measurement in the constraints computation. A complete study of these kind of techniques is developed in [START_REF] Grimson | Object Recognition by Computer: The Role of Geometric Constraints[END_REF], but the complexity remains O((n + 1) m ) in the worst case.

Iterative Closest Point

This algorithm was introduced in [START_REF] Besl | A method for registration of 3-D shapes[END_REF] and [START_REF] Zhang | Iterative point matching for registration of free-form curves and surfaces[END_REF]. It consists of the alternate optimization of the matches and the transformation. Since matches are computed from the distance between points, this algorithm can be easily generalized to features: from an initial transformation f o , one iterates the two following steps:

Matching: Each model feature x i is matched to its nearest neighbor (NN) y j in the scene Y:

π t (j ) = arg min j dist f t x i , y j , i.e., y π t (j ) = NN Y (f t x i ).
Registration: This is usually done using least squares. We have developed in [START_REF] Pennec | L'Incertitude dans les Problèmes de Reconnaissance et de Recalage: Applications en Imagerie Médicale[END_REF] such methods on generic geometric features.

Toward a Generic Framework for Recognition 64

The above matching criterion is generally not symmetric: we can have y j = NN Y (x i ) while the reverse is false. To solve this problem, we have introduced in [START_REF] Pennec | A geometric algorithm to find small but highly similar 3D substructures in proteins[END_REF] a symmetric version of the nearest neighbor where we match x i to y

j = NN Y (x i ) if and only if NN X NN Y (x i ) = x i .
The usual problem of this algorithm is to choose a termination criterion. Since we are using identified features (and not "continuous" features such as points on curves or surfaces), the set of possible matches is discrete, and each interpretation corresponds to a unique transformation. When the algorithm converges, we thus obtain not only the same matches but also the same transformation in successive iterations. This last characterization is computationally cheap and constitutes the termination criterion. (Adding a maximal number of iterations is also useful in case the algorithm does not converge.)

This algorithm is very easy to implement, but it needs a very fast nearest neighbor search in a manifold and an efficient registration computation. Moreover, it can be very sensitive to the initial transformation, and it tries to match the whole model with the whole scene. Thus, unless it is initialized with a very good transformation, it cannot find a very small set of good matches among a lot of clutter (like in proteins). For these reasons, it is often used as a verification algorithm after one of the other algorithms presented here.

Hough Transform

This method was introduced by Paul Hough in 1962 to detect simple geometric objects like lines from points by accumulating evidence in their parameter space. The method was generalized for recognition in the following way: let k be the minimum number of matches needed to compute a unique transformation (or at least a finite number) between the model and the scene. For each k-tuple of matches, we compute this transformation (if it exists), and we accumulate evidence in the transformation space to find the one that maximizes the number of matches (and hence the number of k-tuples matches).

The accumulation stage can be done with a discretization of the transformation space, which serves as an accumulator for the votes cast by the above k-tuple transformations. A sweeping stage is then necessary at the end to find the best transformation hypotheses. One can also keep all the computed transformations in a list and use a clustering algorithm (like in [START_REF] Fischer | An efficient automated computer vision-based technique for detection of three-dimensional structural motifs in proteins[END_REF][START_REF] Pennec | A geometric algorithm to find small but highly similar 3D substructures in proteins[END_REF]). An adaptive subdivision of the transformation space was proposed in [START_REF] Breuel | Geometric Aspects of Visual Object Recognition[END_REF]. At each step, the region is either subdivided or dismissed based on upper and lower bounds on the number of possible matches in the region. However, the linear equations used to efficiently compute these bounds seem difficult to generalize to generic geometric features.

The Hough transform was particularly studied for the inclusion of measurement errors; a synthesis can be found in [START_REF] Grimson | Object Recognition by Computer: The Role of Geometric Constraints[END_REF]. The complexity is O(m k • n k ), but we need to store the transformation space in memory and sweep it at the end of the algorithm to find the maximum. This can be a serious drawback when it comes to 3-D since rigid body motions have a dimension 6 and affine transformations a dimension 9.

Alignment or Prediction-Verification

This technique can be seen as a compromise between interpretation trees and Hough transform: one first hypothesizes a k-tuple of matches that allows one to compute a transformation. Then this hypothesis is verified

VIDERE 1:2
Toward a Generic Framework for Recognition 65 by mapping the model onto the scene and searching for new matches in a given error zone around the model features. The quality score of this hypothesis is simply the number of matches found, or it can be more elaborate. This technique was principally developed in [START_REF] Ayache | Hyper : A new approach for the recognition and positioning of two-dimensional objects[END_REF][START_REF] Huttenlocher | Object recognition using alignment[END_REF][START_REF] Huttenlocher | Recognizing solid objects by alignment with an image[END_REF].

In theory, one should repeat this scheme for all possible k-tuple of matches and keep the best hypotheses. For instance with 3-D points, we need three matches to specify a rigid transformation (k = 3). Thus, there are

C 3 m = m 3 = m!
3!•(m-3)! ways to choose three points among m in the model, C 3 n = n 3 ways in the scene, and 3! to match the two triplets. Therefore, we have m

3 • n 3 • 3! = O(m 3 • n 3
) alignments or hypotheses to verify. In practice, one stops as soon as we estimate that an object is recognized and registered. Moreover, we can use the three invariants of the triplet (the distances between points) to index it in a hash table. This allows one to retrieve in quasi-constant time the compatible triplets in the scene. The complexity falls thus to O(m 3 + n 3 ) for the prediction step. This use of the invariants of the k-tuple can be generalized to geometric features. (See Section 3.6.)

The verification step is particularly important since it should reject the bad hypotheses but keep the good ones. (See Section 6.) It is generally realized with an iterative nearest neighbor with a threshold on the distance for matching. Refining both the transformation and the matches is important if we want to stop as soon as we find a good solution. If we test all hypotheses, we can refine only the best ones.

Geometric Hashing

This technique was introduced by [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]. The idea is to pre-compile the information about model objects in a hash table with an invariant representation (with respect to the action of the transformation group G). This representation should also be redundant and based on local features to allow for the recognition in the presence of occlusions. At recognition time, we need to compute only the corresponding representation of the scene and accumulate evidence for the matching part of the scene with an object.

More precisely, the geometric hashing concerns the case of a subgroup of affine transformations acting on k-D points: one can then define an intrinsic coordinate system (or basis) of a model by choosing at most k + 1 of its points and expressing the coordinates of the others in this basis. For instance, two points are sufficient to define an orthonormal basis in 2-D (in fact a point and a direction are enough), or three nonaligned points in 3-D. We need exactly three points for an affine basis in 2-D and four points in 3-D. The coordinates of other points in this basis are invariant with respect to a global motion of the object.

The idea is to index in a pre-processing step all possible basis in the model by the coordinates of the other model points. (See Figures 4 and5 for examples of hash tables.) At recognition time, one chooses an image basis and computes the coordinates of the other points in it. These coordinates being invariant, we retrieve for each point (thanks to the hash table) the model bases having a point in the same configuration, and we tally a vote for the matching of these model bases with the current scene basis. If the scene basis belongs to an indexed object, the number of votes for the corresponding model basis will be the number of points of the object that are visible in the scene (minus the points of the basis). The complexity is thus O(n k+1 ) in the worst case for the Toward a Generic Framework for Recognition 66 recognition if the access time of a hash table bucket is constant. This is obtained by considering n votes for every of the n k possible scene bases. The verification step is not taken into account here.

The use of a hash table gives a sublinear recognition complexity with respect to a library of objects: one simply indexes all the objects in the same hash table and adds in the information stored to which object each basis belongs.

Geometric Invariant Indexing

The generalization of geometric hashing to geometric features other than points is not so easy, contrary to the alignment or the Hough transform methods. Indeed, the notions of basis and coordinates of other points in such a basis are linked to the structure of vector space, which is generally not available for geometric features. However, one can conceive a similar algorithm using the invariants of a fixed number of features instead of coordinates, the accumulation still taking place in the correspondence space.

If we take for instance 3-D points, the equivalent of the geometric hashing algorithm would be to use in the hash table the invariants of an ordered 4-tuple of points (x i , x j , x k , x l ) instead of the coordinates of x l in a basis constructed from (x i , x j , x k ). For the accumulation, it then seems natural to replace the increment of a basis match by the increment of all four individual matches: if (y i , y j , y k , y l ) is a compatible 4-tuple of points in the scene, we will now vote for each match (x i , y i ).

After this modification, it is no longer necessary to impose the use of a 4-tuple of points: we can use a triplet of points or simply a pair. The complexity is in this case reduced from O(n 4 ) to O(n3 ) or O(n 2 ). However, we must be careful that decreasing the number of points drastically reduces the selectivity of invariants, and the probability of false positives becomes high very fast with the noise. The results of such an algorithm can thus be absurd. (See Section 6.) If we consider a triplet of points, we obtain an algorithm that is very close to the Hough transform, except that the accumulation is in the correspondence space instead of the transformation space. It is possible to combine both approaches by clustering for instance the transformations associated with each possible match. This type of algorithm has been used in [START_REF] Califano | Multidimensional indexing for recognizing visual shapes[END_REF] with complex features (points with local shape descriptors), in [START_REF] Kishon | 3-D curve matching using splines[END_REF][START_REF] Guéziec | Smoothing and matching of 3d space curves[END_REF] for 3-D curves using differential invariants (principal directions and curvatures), in [START_REF] Stein | Structural hashing: Efficient threedimensional object recognition[END_REF] for 3-D range data, and in [START_REF] Fischer | An efficient automated computer vision-based technique for detection of three-dimensional structural motifs in proteins[END_REF][START_REF] Fischer | 3-D substructure matching in protein molecules[END_REF] to reduce the complexity of 3-D pointbased matching from O(n 4 ) to O(n 3 ). In this last case, a pseudo-basis constituted of two points was replacing the standard three points basis for accumulation. We have also used such an algorithm in [START_REF] Guéziec | Medical image registration using geometric hashing[END_REF][START_REF] Pennec | A geometric algorithm to find small but highly similar 3D substructures in proteins[END_REF] with 3-D frame features 3 : we used in this case a pair of frames to index, the invariants of such a pair being the rigid transformation from one frame to another. With this 6-D invariant space, we have a quite selective scheme and a complexity of O(n 2 ) instead of O(n 4 ) for points. However, for the general feature case, the main problem is how to compute the invariants of a k-tuple of features, which will be detailed in Section 5.1.

Error Handling

All these recognition algorithms are quite simple in the theoretical case of exact measurements. In practice, all our measurements are subject to errors due for instance to the deformations of the acquisition system and the image sampling in pixels and gray levels. The algorithms have to be adapted to explicitly take into account the fact that the features are inherently uncertain. For instance, in geometric hashing the measurement error propagates into the computation of invariants, and voting punctually in the hash table (just through the bin corresponding to the computed invariant value) could lead us to miss a great number of matches. We would then obtain false negatives which means that the object is not recognized when it should be. The same thing can happen with the alignment algorithm if the area search for a correspondent is too small with respect to the measurement error of the features enlarged by the uncertainty of the hypothesis transformation.

In this section, we analyze how measurement uncertainty should be taken into account in order to avoid false negatives. This will ensure the correctness of the recognition algorithms: objects will be recognized if they are present in the scene.

Error or Compatibility Zones

To handle uncertainty, we assume that we know an estimation of the measurement errors either as a bound on the values or as a probability law.

Conservative error bounds

The approach of Grimson and Huttenlocher [START_REF] Grimson | On the sensitivity of geometric hashing[END_REF][START_REF] Grimson | A study of affine matching with bounded sensor error[END_REF] or Lamdan and Wolfson [START_REF] Lamdan | On the error analysis of geometric hashing[END_REF] is to propagate an error bound in the invariant computation to obtain an error zone for the votes in geometric hashing or for the search zone during alignment. Using this error zone, we are ensured not to miss a possible match. It is not assumed in this method that the error distribution is uniform but only that the distribution has a compact support that is included in the error zone. To be sure that we do not miss a possible match, we want to compute a conservative error bound (an upper bound) during the operations we perform on features. This fact ensures the correctness of the algorithm. On the other hand, the recursive use of upper bounds ends by giving error zones that are much larger than what is observed.

To illustrate this point, we have computed in [START_REF] Pennec | Correctness and robustness of 3-D rigid matching with bounded sensor error[END_REF] the conservative propagation of error bounds for rigid alignment and geometric hashing of 3-D points (Figure 2). We can see that the predicted error bound is much larger than the values actually observed. In fact, we have established in this case that Mean error = Predicted error 13.2 and Max error < Predicted error 2 .

The conservative bound for this relatively simple example is thus already two times larger than the bound statistically observed. Moreover, the computation of these conservative bounds is especially difficult, and each operation needs a particular hand derivation of the bound.

Probabilistic error Another approach is to consider that we observe the realization of some random vectors and propagate the first moments of the distributions through all computations. This is especially well adapted to introduce a quantitative estimation of the quality of the matches in the matching criterion. This is the core of the probabilistic geometric hashing proposed by Rigoutsos and Hummel [START_REF] Rigoutsos | Massively Parallel Baysian Object Recognition[END_REF][START_REF] Rigoutsos | A Bayesian approach to model matching with geometric hashing[END_REF].

According to the theory developed in the introduction, a random feature is approximated by its mean and its covariance matrix: x ∼ (x, xx ). However, most probabilistic recognition methods assume that the error distribution is Gaussian to combine probabilities in the matching score. This assumption is consistent with our approximation since the Gaussian is the distribution that minimizes the information when we only know the mean and the covariance. 4 However, the main problem of the Gaussian distribution is that its support is infinite. Thus, there is an always non-zero probability of matching any two features. This is catastrophic for the complexity of the matching algorithms since we should theoretically explore the whole correspondence space.

Truncated probabilistic error

The usual solution [START_REF] Rigoutsos | Massively Parallel Baysian Object Recognition[END_REF] is to combine the two approaches by bounding the admissible error using the Mahalanobis distance. (This corresponds to a χ 2 in the Gaussian hypothesis.) The error or compatibility zone around feature x ∼ (x, xx ) is thus:

Z ν (x) = z ∈ M/µ 2 (x, z) = -→ xz T • (-1) xx • -→ xz ≤ ν 2 , ( 4
)
where ν 2 is a (generally global) threshold that can be interpreted in the Gaussian case as a χ 2 . It is interesting to compare with the previous bounded error model:

Z ε (x) = z ∈ M/ dist(x, z) 2 = -→ xz T • -→ xz ≤ ε 2 (5)
where the threshold ε is metric and thus harder to choose than a threshold without dimension. The truncated probabilistic model can thus be seen as a generalization of the bounded error model where the information matrix (-1) xx is used as a local metric. Moreover, one can show [10, chap. 5, p. 152] that propagating covariance matrices using the Jacobians can be interpreted in a deterministic way as the first-order approximation of the error bound propagation. [START_REF] Boissonnat | VIDERE 1:2 Toward a Generic Framework for Recognition 85[END_REF]. By the way, it is using this very property that we have defined the Gaussian distribution on a manifold in [START_REF] Pennec | L'Incertitude dans les Problèmes de Reconnaissance et de Recalage: Applications en Imagerie Médicale[END_REF].

VIDERE 1:2

Toward a Generic Framework for Recognition 69

Handling Uncertainty in Matching Algorithms

We describe in this section how to modify the recognition algorithms in order to make them handle an a priori known uncertainty on features. It is, however, desirable to verify a posteriori after recognition and registration that our guess about feature uncertainty is good. One can also think of iterating the whole process (matching, registration, and estimation of the noise on features).

Interpretation trees

To handle uncertainty in this algorithm, we just have to propagate the uncertainty in the computation of unary, binary, and higher-order invariants and perform the constraint satisfaction test using a Mahalanobis distance between uncertain invariants.

Hough transform Assume that we have found a k-tuple of model features and a k-tuple of scene features that have compatible invariants. (This can be done using an invariant indexing technique.) We compute the transformation f between these k-tuples and its uncertainty ff . We obtain in fact a probabilistic transformation f ∼ ( f, ff ). To verify if these k matches are correct, we use the Mahalanobis distance between matched features after registration; but since the transformation is computed by minimizing these distances, we should not take into account the uncertainty of the transformation to obtain an unbiased result. (Otherwise the Mahalanobis distances would be underestimated.) This means that we test if µ 2 ( f x i , y i ) ≤ ν 2 for each of the k matches.

Since we need at least k matches to compute the transformation, this one is not reliable if at least one of the k tests fails. We reject in this case the transformation. If all tests are passed, the transformation is a possible one, and it is added in the accumulator with its uncertainty.

There are two main techniques for the accumulation. The most widely known is to vote for the bin where our transformation falls in the sampled transformation space. This is the origin of the term "accumulator." With the error handling, we now have to vote for all the bins of the accumulators that intersect the compatibility zone Z ν (f ) of our probabilistic transformation. This problem is not so easy to solve efficiently and will be detailed in Section 5.4. We also have to sweep the sampled transformation space at the end of the algorithm to find the maximum score transformations. An alternative is to maintain during the algorithm a list of the best transformations. (The question is then how many do we have to keep.)

The second technique is to add each possible transformation f to a list and afterward use a clustering algorithm to extract and merge the sets of compatible transformations. We used a similar technique in [START_REF] Pennec | A geometric algorithm to find small but highly similar 3D substructures in proteins[END_REF].

Alignment

The computation of the transformation between k features and the verification of theses k matches is exactly the same as for the Hough transform. The difference is that we now verify the hypothesized transformation by looking for additional matches with the other features. In this process, the uncertainty of the transformation has to be taken into account. One then searches for each model feature x i the nearest neighbor y j = NN(f x i ) in the scene (according to the canonical or the Mahalanobis distance and possibly with the symmetry constraint on the nearest neighbor). The match is accepted if the following test is passed:

VIDERE 1:2
Toward a Generic Framework for Recognition 70

µ 2 (f x i , y j ) < ν 2 .
If the number of matches is sufficient for the hypothesized transformation, it is interesting to verify the global consistency of matches and an improvement of the transformation. We can indeed seriously reduce the transformation uncertainty using all matches.

ICP and global consistency verification

If we have the matches (as with alignment), the following step is to recompute the transformation minimizing the Mahalanobis distance between all these matches. To improve robustness, we can rule out the outliers by verifying once again the Mahalanobis distance of the matches, but, since the new transformation is computed from these very matches, one has to do it without the transformation uncertainty to be unbiased. The test is thus for each match:

µ 2 ( f x i , y j ) < ν 2 .
The process can be iterated until convergence. To realize a real ICP from this verification scheme, it is sufficient to recompute the nearest neighbor of each model feature, keeping only those who pass the above test.

Hashing and geometric invariant indexing

For geometric invariant indexing, we need to propagate the feature uncertainty in the invariant computation (as for the interpretation trees), and index these probabilistic invariants with their error zone. This last step raises the same problems as for the indexing of probabilistic transformations in the Hough transform. It will be detailed in Section 5.4.

Handling error in the geometric hashing algorithm is more complex. Consider for instance the case of 3-D points: we need first to fix the way we compute the basis from three points and estimate the uncertainty of this local basis (i.e., the uncertainty of the transformation f from the canonical basis to the local one) with respect to the uncertainty of the three points. The best way to do this is unclear, even if a least squares (between the three points and three fixed points in the canonical basis) seems to be adapted. We have in this case already developed the algorithm to get the registration and its uncertainty. Computing the coordinates of other points in the local basis then corresponds to the action of the probabilistic transformation f . In this process, we must be careful to take into account the uncertainty of the transformation while computing the uncertainty of the invariant coordinates (otherwise a slightly large error of the basis features would rule out most of the matches). The last modification is to index (and retrieve) as above our probabilistic invariant coordinates with their error zones.

Some Key Problems

Up to now, we have formalized the main matching algorithms used on points in terms of features and seen how to modify them in order to explicitly handle uncertain measurements. In this section, we investigate some key problems that are raised by these algorithms and that constitute the difficult part to implement either due to the fact that we are using features or simply because of the uncertainty.

We have identified four main problems. For ICP, we know how to compute the transformation and its uncertainty, but the problem is to find efficiently the closest neighbor in a manifold (with a non-Euclidean Riemannian metric). A related problem for the Hough transform is to Toward a Generic Framework for Recognition 71 cluster the transformations. Of course, the clustering of geometric features can be of use for many other important algorithms. The other algorithms (including the Hough transform) use more or less unary, binary, or higher-order invariants. The problem is not only to compute them but also to compare them efficiently and reliably. Last but not least, hashing is widely used to improve the algorithmic search for compatible features or invariants. How can we do this in presence of error and in a manifold?

n-ary Invariants: Shape Space

We have previously talked about unary, binary, or higher-order invariants as the characteristic invariants of the shape of k ordered features. If we can easily see that this notion corresponds to the distance for a pair of (rigid) points and, for instance, the three interpoint distances in a triplet of (still rigid) points, it is much harder to imagine what are the similarity invariants of a 4-tuple of points or the affine invariants of five oriented points.

We believe that an approach based on the shape theory is well adapted to tackle this problem. This theory was principally developed by Kendall and Le [START_REF] Kendall | A survey of the statistical theory of shape (with discussion)[END_REF][START_REF] Le | The Riemannian structure of Euclidean shape space: A novel environment for statistics[END_REF] on points under similarities and rigid transformations. The idea is to characterize the configuration space of a set of k ordered features. By configuration, we mean what is invariant: the shape. The method is the following: a k-tuple of features is an element of M k . To obtain its shape, we identify all k-tuples that can be identified using a suitable transformation f ∈ G. The k-shape space is thus the quotient space I k = M k /G that Kendall denotes by (M, G, k). One can also see this as the "factorization" of a k-tuple of features into a shape i ∈ I k and a "position" f ∈ F ⊂ G of this shape in space.

With this approach, the first problem is to find a metric on the shape space compatible with our metric on the manifold M and on the group G. Roughly, this means that we want to have the same measurements on invariants independently of the position of the original k-tuple in space, and this is equivalent to finding a "factorization" of the metrics such that M k = G × I k . With the metric on the invariant space, we can determine geodesic and the exponential chart (the "most linear chart" with respect to the metric) at each point of the manifold.

Then, from a theoretical point of view, we can compute the mean invariant and its covariance matrix and the Mahalanobis distance. From a computational point of view, we need to find a way to implement the exponential chart at each point of I k as we no longer have a placement function as in the homogeneous manifolds to identify things at feature x with things at the origin. The problem is however simpler than for features, since the basic operations on probabilistic invariants are reduced to Translation between a k-tuple of features and the pair k-shape/position:

X ∈ M k ↔ (i, f) ∈ I k × G. Distance between two k-shapes: dist(i 1 , i 2 ).

Mahalanobis distance between two probabilistic k-shapes: µ(i 1 , i 2 ).

From these basic operations, we can construct algorithms to merge invariants, estimate the mean, etc. For efficient matching algorithms, we also need to find the nearest neighbor in this manifold or index these uncertain invariants.

Clustering of Uncertain Geometric Features

We are interested here in the clustering of transformations for the Hough transform. However, since the transformation group is a manifold, a more general problem is the clustering of geometric features on a manifold. The problem is somehow simplified with respect to the standard clustering problem since we have an uncertainty estimation on our features: we know thus the "scale" of the clusters. On the other hand, we generally do not know the number of classes.

It is possible to generalize some classical techniques based on the distance of points (see for instance [START_REF] Duda | Pattern Classification and Scene Analysis[END_REF][START_REF] Jain | Algorithms for clustering data[END_REF]) by replacing the distance between points with the Mahalanobis distance between uncertain features. On the contrary, space sampling techniques are less likely generalizable since we no longer have a vector space. (See also Section 5.4.)

We have used in [START_REF] Pennec | A geometric algorithm to find small but highly similar 3D substructures in proteins[END_REF] a very simple but rough technique: the information (or its opposite value, the entropy) of a random feature f is related to the log of the determinant of its covariance matrix. We can thus choose the most informative feature among the set to cluster and iteratively merge the closest feature to the current state estimate (according to the Mahalanobis distance). Each used feature is removed from the set. Once there are no more features to merge, we have obtained one cluster represented by its mean feature and we iterate the clustering stage on the remaining features. In this process, an efficient way of finding the nearest (Mahalanobis) neighbor would be an important improvement for the complexity. A more rigorous algorithm would be to let the different clusters compete with each other and to compute the mean feature of a cluster with a Mahalanobis distance minimization at each step.

Nearest Neighbor in a Manifold

This is an old problem and it is well studied on points in computational geometry. The corresponding notion is the Voronoï diagram. However, constructing and using this diagram is complex, so most techniques rely on space-partitioning methods such as k-D trees [START_REF] Preparata | Computational Geometry, an Introduction[END_REF].

This last technique iteratively subdivides the space along each axis in turn and relies on the equivalence between the Euclidean L 2 norm and the L ∞ norm (maximum coordinate) which is separable along the axes. This allows us to give upper and lower bounds on the Euclidean norm with respect to the maximal difference between coordinates. Its generalization to a Riemannian manifold seems difficult since there no longer is a global coordinate system to define a L ∞ norm and alternate the search along the axes.

The Voronoï diagram can be generalized to a Riemannian manifold, but very little work exists. One can cite [4, chap. 18] for a hyperbolic manifold, but the determination of the diagram relies on the fact that the manifold has a negative curvature and thus there exists a global diffeomorphism with R n (in this case a projection). These kinds of techniques are not applicable to positively curved manifolds such as spheres or projective spaces (which includes 3-D rotations). However, other techniques are possible. Watson developed a method to compute the Voronoï diagram on spheres S n with their canonical Riemannian metric. Using the very special properties of the exponential charts, we think that it is possible to compute the Voronoï diagram on any homogeneous manifold with an invariant Riemannian metric. However, the efficiency of such a construction for the nearest neighbor problem is not ensured. To conclude on this problem, we note that none of these techniques can apply for the nearest neighbor according to the Mahalanobis distance. Indeed, they rely on data preprocessing that makes use of the metric, which is only known at query time with the Mahalanobis distance. (It uses the covariance of the query feature.)

Uncertain Indexing

In all the algorithms that make use of invariants we have the same problem: how to retrieve efficiently the nearest or compatible invariants among a predefined set. The "brute force" method is to compare our query invariant with all the others, which has a linear complexity with respect to the number N of indexed invariants. There is no preprocessing stage, and the memory requirement is also O(N ). There are basically three techniques to reduce the complexity. We have investigated in Section 5.3 the Voronoï diagram to find the nearest neighbor. We can subdivide the space sampling methods in two: either the sampling depends on the data, as in k-D trees, which generally gives a search time of O(log N), or the space is sampled in a fixed way and the contents of a bin can be retrieved in quasi-constant time with a hash table . In this section we focus on this last technique.

In an ideal world with no noise on the data, we just have to determine a space sampling and a hash function that transform the coordinates of a bin into a code. This code indexes (through an array) the list of data having the same code. This hash function should scatter as much as possible the codes to obtain few empty lists and a very short mean list length. (Ideally there should be one and only one data per code.)

The introduction of the error on our invariants raises some problems in this algorithm: we now have to index or retrieve them using an error zone that generally intersects several bins in the sampled space. We have first to determine which bins. Next, all the bins create supplementary entries in the hash table. We will no longer have N entries but much more, for instance around 6N for the 2-D hash table of Figure 4. It is then important to quantify the mean number of indexed bins per data to have an idea of the mean number of objects by index (the mean length of the list of objects for a given code). The real complexity of hashing is proportional to these two factors. Since the hash function can associate the same code to very different locations and the volume of the compatibility zone is smaller than the indexed volume, it is moreover necessary to verify the compatibility of the query feature and those found in the corresponding bins of the hash table.

In [START_REF] Pennec | Correctness and robustness of 3-D rigid matching with bounded sensor error[END_REF], we have analyzed the complexity of the different types of hashing used by recognition algorithms on points. Let us index for instance points of R d with an error zone bounded in norm by ε and a Cartesian sampling with width l (such as in Figure 5). The mean number of indexed bins by point (i.e., intersecting the error zone) is about n = (1 + 2 • ε/ l) d . Using this value in a false positives analysis, we concluded that l ε was a good trade-off. Thus, for 3-D points, we have a mean number of nine indexed bins. For a fixed dimension, this is only a constant multiplicative factor in the complexity of the algorithm, but this factor is exponential with respect to the dimension. Hence, for the (rigid) geometric hashing of 3-D points, using the three invariants of the basis (the distance between the three points) and the coordinates of the fourth point in this basis, we index in a 6-D space and we obtain a mean number of index bins per point of 729! Moreover, even if we assume the same metric error bound on all points, its propagation through the computation of invariants gives a bound that depends on the invariant position, as can be seen in Figure 5. Therefore, an adaptative space sampling is sometimes preferable.

Up to now, we have only raised the problems for the simple invariants of points. If we consider now a frame pair, the associated binary

Toward a Generic Framework for Recognition 75

invariant is also a frame. (It is in fact the rigid transformation from one frame to the other expressed in one of the frames.) The invariant space is this time SO 3 × R 3 and no longer R d . The first question is how to sample regularly a space closing on itself such as the sphere or the set of rotations SO 3 . This problem is linked to the way we compute the intersection of the error zone with the bins, which has to be very efficient for the hashing to remain interesting.

The adaptation of hashing techniques to uncertain geometric features and invariants is thus a difficult problem, and, unless new techniques were to be designed, the speedup they were used for is not always conserved. Efficient searching for compatible features or invariants is, however, a crucial problem for the efficiency of recognition algorithms based on geometric features.

Performance Analysis: False Positives

With the correct handling of measurement errors, we have ensured the correctness of our recognition algorithms: there will be no (or very few) false negatives. On the other hand, we now have a larger probability of false positives (or phantom recognition). Indeed, using an error zone instead of a unique feature allows us to match features that fall by chance in this area, and, when this probability is sufficiently high, individual false matches can combine themselves (conspiracy) to produce an important matching score. We then will estimate that we have recognized an object that is not present.

There are two principal sources of false positives. The first one is that we are looking for a local consistency that can be insufficiently constrained. For instance, this is the case when we use only unary and binary invariants to predict global matches. A verification step is thus needed to insure the global consistency of the matches. This verification step can moreover improve the registration accuracy and rule out outliers. Since this step can be computationally expensive, we should limit the number of false positives as much as possible.

The second source of false positives is less obvious and comes from the simple modeling of the image (or the data) by features. Some important information can be dismissed in this process, and we can end up with a correct recognition from a feature point of view (a globally consistent matching) that is incorrect from the image or data point of view. This problem is inherent to the ascendent organization of information in image processing but can be minimized by using more-adapted and more-informative features. For instance, we will see in Section 6.2 that adding trihedra to 3-D points to form frames allows us to get rid of more than 80% of the individual false matches even with a very noisy trihedron (a standard deviation of 90 • !). In real applications (registration of medical images from "extremal points" [START_REF] Pennec | A framework for uncertainty and validation of 3D registration methods based on points and frames[END_REF]), the standard deviation of the trihedron is about 10 • , and the probability of false positives using frames drops off drastically with respect to points (Figure 7). In other experiments in molecular biology [START_REF] Pennec | A geometric algorithm to find small but highly similar 3D substructures in proteins[END_REF], we observed that using frames instead of points could rule out some biologically nonsensible matches.

In the following, we present a new method to analyze the probability of false positives for some matching methods. Then, we tackle a more general problem: the computation of the intrinsic complexity of the recognition problem, independently of the method used. Toward a Generic Framework for Recognition 76

Problem Formulation

The probability of obtaining a false positive has been mainly studied in [START_REF] Grimson | On the sensitivity of the Hough transform for object recognition[END_REF][START_REF] Grimson | On the verification of hypothesized matches in model-based recognition[END_REF][START_REF] Lamdan | On the error analysis of geometric hashing[END_REF][START_REF] Grimson | A study of affine matching with bounded sensor error[END_REF] for recognition problems from a library of exact models based on points and lines. Here, we still consider that the model is exact but based on m geometric feature instead of points. The scene consists of n noisy features. In the model and in the scene, τ features are in the same configuration (up the noise). The (mτ ) and (nτ ) other features are supposed to be randomly distributed in the model image I m and the scene image I s . In fact, if the meaning of this sentence is quite clear for points, it needs to be detailed for features. By "image," we mean the set of possible measurements of our features in the real image. For instance in a 3-D image, the position of a frame is constrained to be in the image volume U ⊂ R 3 but the trihedron part could be anything. Thus, the image (from a frame point of view) is I = SO 3 × U ⊂ M. The random distribution of features in an "image" I is of course the uniform (i.e., invariant) distribution on this set.

To compute the false positives probability, we also need some hypotheses about the matching criterion. Here, we use one of the simplest: the criterion is the number of matches. A set of matches is accepted if its score after verification is greater than a given threshold. The matching algorithms being modified in the previous section not to miss a match, we can take τ for this threshold.

Selectivity: Probability of a (Single) False Match

Let f be a hypothesized transformation from the model to the scene: the exact feature x is matched with the random one y if the transformed feature f x is in the error zone Z(y). This error zone can be based on a truncated probabilistic model (Equation 4) or on a bounded error model (Equation 5).

We call selectivity the probability of accepting this match if one of the features is an outlier (its position is random and uniform in the image). By symmetry, we assume that feature x is the outlier. Then, the probability of a false match is the conditional probability:

P (f x ↔ y) = P (f x) ∈ Z(y)|x ∈ I m .
The uniform distribution on the set I m is given by the invariant measure dM on the feature manifold (see [START_REF] Pennec | Uniform distribution, distance and expectation problems for geometric features processing[END_REF] for a discussion), normalized by the volume of the set V(I m ) = I m dM. Thus, the selectivity is

P (f x ↔ y) = (f I m )∩Z(y) dM V(I m ) = V (f I m ) ∩ Z(y) V(I m ) .
If the volume V(Z(y)) of the error zone is small with respect to the image volume V(I m ), we can consider that the transformed image f(I m ) either contains the whole error zone or does not intersect it at all. This allows us to approximate the above probability by

P (f x ↔ y) = ε V Z(y) V(I m ) where ε = 1 if f (-1) y ∈ I m ε = 0 otherwise.
A desirable property for our "error volume" Z(y) is to be comparable at every point since we usually fix the same bound for error on all the points. This means that, for any feature y , there exists a transformation f such that y = f y and Z(y ) = f Z(y). The error volume is

VIDERE 1:2
Toward a Generic Framework for Recognition 77 said to be homogeneous. A stronger hypothesis is that for every transformation f, the error volume on the transformed point is the transformation of the error volume: Z(f y) = f Z(y). The volume is said to be isotropic in this case, and is completely determined by its shape around the origin. (See [START_REF] Pennec | L'Incertitude dans les Problèmes de Reconnaissance et de Recalage: Applications en Imagerie Médicale[END_REF] for an analysis of noise models.) In both cases (homogeneity and isotropy), the volume of the error volume is invariant and can be computed at the origin. In the case of a homogeneous probabilistic model, it depends only on the covariance at the origin = J (f y ) (-1) • yy • J (f y ) -T :

V Z ν (y) = V 0 = x T • • x≤ν 2 dM( x).
Example with frames To keep this example as simple as possible, we consider here a bounded error model on the position and a separated one on the orientation. Thus, two frames are matched if the distance between their points is less than a threshold d 0 and if the rotation needed to adjust their trihedra has an angle less than a threshold θ 0 (this angle is θ = r (-1) x

• r y ). This error zone is isotropic. The volume is thus invariant and we can compute it at the origin: a frame f = (r, t) is in the error volume Z(I d) if θ = r < θ 0 and t < d 0 . Using the invariant measure on rigid transformations dM(r, t) = sin 2 ( r /2) r 2 dr dt, we can compute the volume of the error zone:

V 0 = θ<θ 0 t <d 0 dM(r, t) = θ<θ 0 sin 2 (θ/2) θ 2 dr • t <d 0 dt V 0 = 2π(θ 0 -sin(θ 0 ) • 4π 3 d 3 0 .
If we assume a cubic image of side l (256 for instance), this gives a Euclidean volume V I = l 3 for points in which trihedra are not constrained: the rotation volume is 2π 2 . Finally, we obtain the basic probability of false match:

P (f x ↔ y) = ε • η with η = θ 0 -sin θ 0 π 4 3 d 0 l 3 .
We have isolated in the first term the probability of false match due to the trihedra only, which reflects the gain in selectivity when using frames instead of points. This function is plotted in Figure 6 and shows very interesting results: even for a bound of θ 0 = π/2 = 90 deg, more than 80% of the random matches are rejected. For a more realistic bound of θ 0 = π/10 = 18 deg, the probability of a false match drops to 0.0016: we would have to divide the bound on the position by 10 to obtain an equivalent selectivity using points only.

We will see in the rest of this section that it is sometimes useful to assume that the image is spherical, for instance with a diameter d = √ 3 • l. The volume of the (frame) image in this case is V = 4 • π • (d/2) 3 /3 and the selectivity becomes:

η = θ 0 -sin θ 0 π 2 • d 0 d 3 . VIDERE 1:2
Toward a Generic Framework for Recognition 78 Figure 6. Basic probability of a false match for trihedra with a bound on the angle for the adjustment rotation of theta. Since the formula of the selectivity is multiplicative for frames, this curve is also the gain in selectivity when using frames instead of just points (the ratio of the selectivity of frames over the one of points). theta In this section, we investigate some basic false positives analysis for the Hough transform and the alignment algorithm. Since many variants of these algorithms exist, the purpose is not to give a precise number of false positives but to explain with practical examples what techniques can be used. These techniques will be generalized in the next section to the analysis of the intrinsic number of false positives, independently of the algorithm used.

Probability of Choosing a "Correct Basis" in the Model

For the Hough transform and the alignment method, we begin by finding the possible matches of k features (k being the minimal number of matches to determine uniquely a transformation). By reference to the alignment algorithm, we call it a hypothesis. A hypothesis is correct if the k matches are all correct. One false match is sufficient to give an incorrect transformation.

There are A m k = m k • k! possible sets of k ordinated features in the model, but only A τ k are correct. Thus, the probability of choosing an incorrect "basis" is (with an approximation for k τ ):

p (m,k) = 1 - A τ k A m k = 1 - τ ! • (m -k)! m! • (τ -k!) 1 - τ m k .

Number of "Compatible Bases" in the Scene

Assume that we have chosen a basis in the model. A scene basis is compatible if there exists a transformation f that superimposes the k model basis features with the k scene basis features.

Let us investigate the case of a fixed (or known) transformation f. Assuming that the distribution of scene features is uniform, the probability that no features fall in one of the error zones of the transformed model basis is (1ε • η) n . Thus, the probability of finding at least one match in the scene for each model basis feature (under the transformation f) is

p k (f) = k i=1 (1 -(1 -ε i • η) n ) with ε i = 1 if f x i ∈ I s 0 otherwise.
Now, we just have to integrate over all possible transformations to obtain the mean number of compatible bases in the scene. We note that VIDERE 1:2 Toward a Generic Framework for Recognition 79 the expression above is null if one of the ε i is null. Thus, we can rewrite it as

p k (f) = (1 -(1 -η) n ) k
i=1 ε i and the integral becomes:

P k = f∈G p k (f) • dG = (1 -(1 -η) n ) k • f∈G k i=1 ε i (f) • dG.
An upper bound of the value of the right term α is obtained by only asking for the intersection of the images:

k i=1 ε i (f) = 1 ⇒ f I m ∩ I s = ∅.
In the case of 3-D rigid transformations, we can obtain a rough estimation of this bound as follows: let d be the image diameter (for instance d = √ 3 • l for a cubic image of side l). Assuming the images are spherical, and taking the origin at the center, we can make any rotation followed by any translation of length less than d and the images still intersect. On the other hand, any translation of length greater than d separates the two images. Thus, we can bound our integral by

α ≤ 2 • π 2 • t <d dt = 8 3 • π 3 • d 3 = (2 • π • d) 3 3 = (2 • π • l) 3 .
Thus, given a model basis, the mean number (nonnormalized probability) of corresponding bases in the scene image is 3 3 .

P k = α • (1 -(1 -η) n ) k = α • (n • η) k with α ≤ (2 • π • d)

Hough Transform: Proportion of False Transformations

Ideally, we should compute the mean number of true and false hypotheses in each bin of the transformation space. To simplify the problem, we just compute the percentage of good hypotheses. Roughly, this amounts to assuming a uniform distribution of both good and false hypotheses in the transformation space.

We have A τ k correct hypotheses. But for each (correct or incorrect) model base, we find a mean number of P k false hypotheses, which give P k • A m k false hypotheses. Thus, the proportion of false transformations in the accumulator is

P H ough = P k • A m k A τ k + P k • A m k = 1 + A τ k P k • A m k (-1) 1 + 1 α τ η • m • n k (-1)
.

Probability of Accepting a Verification

For the alignment algorithm, we have to verify each hypothesis. Given a hypothesis, k features are already used to constitute the scene basis. Thus, there are nk error zones in the scene. The probability that one of the mk model features avoids all the error zones is (1η) n-k . Thus, the probability that it falls in at least one the error zone is

p = 1 -(1 -η) n-k (n -k) • η.

Toward a Generic Framework for Recognition 80

The probability to match exactly j features among the mk model features is the binomial

B (m-k,p) (j ) = m -k j p j (1 -p) m-k-j .
Thus, the probability of accepting a hypothesis is the sum of the probabilities to have more than τ matches:

q = j ≥τ B (m-k,p) (j ) = 1 - τ j =0 m -k j p j (1 -p) m-k-j = I p (τ + 1, m -τ ),
where I p (a, b) is the normalized incomplete Beta function.

Assuming that p 1, we can approximate the binomial B (m-k,p) by the Poisson distribution of parameter λ = (mk) • p:

q 1 -e -λ τ j =0 λ j j ! with λ (n -k) • (m -k) • η.

Alignment: Proportion of False Positives

We have A m k possible bases in the model, and a mean number of P k false hypotheses for each one, each hypothesis being accepted with a probability q. Thus, we accept a mean number of A m k • P k • q incorrect matchings. We also have A τ k correct hypotheses that will give correct matchings. Thus, the proportion of false positives is

P align = q • P k • A m k A τ k + q • P k • A m k = 1 + A τ k q • P k • A m k (-1) 1 + 1 α • q τ η • m • n (-1)
.

The probability q can be directly interpreted as the filtering ratio of the verification step.

Intrinsic Probability of False Positives

In the previous section, we try to characterize the way two algorithms act by propagating step by step the probability of incorrect matches. One can also ask for the intrinsic probability of false positives, independently of the algorithm used: what is the probability to obtain a score τ if the model and the scene are in fact constituted of m and n features randomly (uniformly) distributed in the images?

In fact, we have developed in the previous section the basic tools to answer this question. Indeed, finding τ matches means that there exists a transformation f such that τ model features fall in the error zones of the scene features. Thus, we just have to compute the probability of obtaining τ matches for a given transformation, and then integrate this probability for all possible transformations.

We have n error zones in the scene. Given a transformation f, the probability for one model feature to fall in at least one of these zones is p η • n. Thus, the probability that exactly j of the m model features fall Qualitative estimation of the number of false positives involving at least x matches in MR images of 2,500 features. Comparison between frames and points: we need roughly five times more point matches than frame matches to obtain the same probability (10 frames and 56 point matches for a probability of 10 -10 ). The actual match found involves about 500 features, and its probability of being a false positive is thus practically zero. in a scene error zone is the binomial B (m,p) (j ). Therefore, the probability of matching at least τ of them is

q = j ≥τ B (m,p) (j ) = I p (τ + 1, m) 1 -e -λ τ j =0 λ j j ! , with λ = m • p = n • m • η.
To obtain the mean number of false positives, we just have to integrate over all possible transformations, which amounts to multiplying it by α. Finally, with our spherical image of diameter d, the probability of matching τ features among m and n random features is

= α • q (2 • π • d) 3 3 •   1 -e -λ τ j =0 λ j j !   with λ = n • m • η. ( 6
)
This equation is valid only when λ = n • m • η 1. Moreover, as our analysis is very qualitative, we can trust it only for very small values. However, notice that the influence of the integration over admissible transformations on the number of false positives is linear. This is to be compared with the exponential influence of the selectivity η in the accepting probability q. Hence, the estimation of α does not need to be really accurate (or could be adjusted a posteriori from experiments). This effect is visible in Figure 7: the difference in selectivity between points and frames induces very different false positives probabilities, whereas a variation on α would produce only a small vertical shift of the curves in logarithmic scale.

Application Example in Medical Imaging

In the case of the MR images of the associated demo, the volume images are 256 • 256 • 165 mm and we typically extract 2,500 extremal points (modeled by semi-oriented frames). Among them, about 500 are matched in a registration.

Frame selectivity Combining the error on the model and on the scene features, we obtain a roughly diagonal covariance matrix: On such a small error zone, we can approximate the invariant measure on frames by the Lebesgue measure:

sin 2 (θ/2) θ 2 • dθ • dt (1 + O(θ 2 )) • dθ • dt 4 .
Thus, the volume of the error zone defined by a χ 2 limit of ν 2 is

V 0 = x T • (-1) •x≤ν 2 dx 4 .
To compute this integral, we make the change of variable x = ν • • y, where is a square root of . We get

V 0 = 1 4 • y T •y<1 ν 6 • | det( )| • dy = ν 6 • det( ) • π 3 6 • 1 4 .
For a χ 2 bound of ν 2 = 16, this gives V 0 1.56. The total volume of the image is 256 • 256 • 165 = 1.0810 7 for the frame position and r ≤π sin 2 ( r /2) r 2

• dr = 2 • π 2 for the frame trihedron. In fact, extremal points are modeled using semi-oriented frames, and we have to divide this volume by 2. Thus the image volume is

V(I) = 256 • 256 • 165 • π 2 = 1.067 10 8
Thus, the selectivity is:

η = V 0 V(I) = 1.46 × 10 -8 .
Point selectivity In the case of points, the noise model becomes isotropic with a standard deviation σ = 0.50. We use a 3-D χ 2 of ν 2 = 6. The volume of the error zone is now

V 0 = ν 3 • σ 3 • 4 3 • π 21.7
, and the image volume is V(I) = 1.08 × 10 7 . Thus, we obtain a selectivity of η pt = 2.01 × 10 -6 . This is a loss of two orders of magnitude with respect to frames!

Probability of false positives

To obtain an upper bound on the number of false positives, we take as image diameter d = √ 3 • l with l = 256, which gives an integration factor α = (2 • π • l) 3 = 4.16 × 10 9 . The number of features is n = m = 2,500. We plot in Figure 7 the probability of a random match of τ features using Equation 6. In both cases (frames and points), the probability of a false positive involving 500 matches (the observed value in real images) is null at the machine accuracy. However, the difference between points and frames is visible for smaller sets of matches. For instance, if we decide that a match is not a false positive if its probability to be so is less than = 10 -10 , we accept matches made of at least 10 frames or 56 points. In fact, we need roughly five times more matches using points than using frames to obtain the same probability of false positives! The interest of using frames instead of points is clear.

Moreover, complex geometric features do have more invariants than simple features, which can considerably reduce the complexity of matching. We were for instance able to reduce the complexity of the 3-D substructure matching problem from O(n 4 ) to O(n 2 ) using frames instead of points. Using geometric features more informative than points is thus crucial to increase the robustness of recognition algorithms and reduce their complexity.

Application Example in Molecular Biology

In [START_REF] Pennec | A geometric algorithm to find small but highly similar 3D substructures in proteins[END_REF], we use frames to model amino acids in proteins. (See the html demo.) In the example of this article, the radius of the protein is about 20 Å, which gives an "image volume" of V I = 4 3 • π • r 3 • (2 • π 2 ) 661500 for frames and V I = 4 3 • π • r 3 33510 if we use just points. We use a diagonal covariance matrix with standard deviations σ r = 15 o = 0.15 rad for the trihedron and σ t = 0.35 Å for the position. The χ 2 values are 16 for frames and 8 for points. Thus, the volume of the error zones are V 0 = χ 6 • σ The integration factor alpha is α = (8 • π • r) 3 /3 = 4.2 × 10 7 , and we have 65 and 105 amino acids in the two proteins. With these values, we need to match at least 5 frames or 20 points to obtain a probability of false match less than 10 -10 , and the main substructures we found in common score 22, 16, and 13 matches. In this case, using frames is critical to ensure the significance of our matching results.

Discussion: Using the False Positives Analysis to Improve the Algorithms

When we tackle a new problem, the first application of this analysis is to provide the intrinsic complexity of the recognition problem. Depending on the type, the error model, and the number of features used, we can compute the threshold τ on the number of matches to use in order to have a given probability of false positives (for instance = 10 -10 ). If this threshold is well above the number of matches we expect, we can proceed and choose an algorithm with a low-complexity basis. Otherwise (as in the above example with proteins), we expect any algorithm to reach its maximal complexity and possibly give false positives. In this case, a solution is to model the images by more-complex geometric features, i.e., to take into account more information in the feature-based recognition process. However, we should notice that the false positive analysis we presented is based on several limitative assumptions (exact model, uniform distribution of features, etc.) which are hardly ever verified in real cases. For instance, in medical images of the head, extremal points are not uniformly distributed in the image, but more or less uniformly distributed on the surface of the brain and the skull, which can be roughly considered as spheres. Thus, the performance analysis we presented should not be considered as quantitative estimations but only as indicative of a general behavior.

A very interesting extension would be to compute the probability of false positives online, during the recognition algorithm itself. This would allow us to take into account the specific distribution of the model and scene features. In our generic framework for feature-based recognition, this could be done by adding a basic operation on features computing the selectivity of matches using the covariance matrix and other operations to propagate this probability through the computations, as we did for the uncertainty. Doing so, we could evaluate the quality of a match at any stage of the algorithm and give up or delay the processing of bad hypotheses. This should be related with probabilistic matching techniques developed in [START_REF] Rigoutsos | A Bayesian approach to model matching with geometric hashing[END_REF].

Toward a Generic Framework for Recognition 84

As far as a specific recognition algorithm is concerned (for instance the Hough transform or the alignment method), we could extend the false positive analysis to compute the mean number of hypotheses to verify before finding a good one, the mean complexity of the algorithms, etc. However, these computations rely on the specific settings of these algorithms, and we cannot develop the analysis for all approaches: we leave it to the reader to analyze his own algorithm. We believe that such an analysis may help to discriminate between different options or settings of an algorithm depending on the application context.

Conclusion

We have formalized in this article the main matching algorithms in terms of geometric features and showed how to modify them in order to incorporate explicitly and rigorously the uncertainty of measurements. The drawback of the correctness of the algorithms is the presence of false positives. We developed a new method to analyze the probability of false positives in two particular algorithms and generalized it to the evaluation of the intrinsic complexity of the matching problem, independently of the method used. Doing so, we showed that using more-informative features, such as frames instead of points, can drastically reduce this probability and allow for the use of lower-complexity matching methods.

For a computational point of view, we have identified four still-open problems to implementing these generic matching algorithms. The first one is linked to the structure of the space of n-ary invariants, and we have formalized it with the shape-space theory. We believe that it is not only possible to characterize, in that way, the manifold of invariants but also a suitable metric structure. The second problem concerns the clustering of features on a manifold, and more particularly with uncertainty information. The method we use to solve this problem can certainly be improved. The last two key points are linked to the efficiency of searching algorithms. This is for the first part the search of the nearest neighbor, with respect to the Riemannian metric or the Mahalanobis distance, and for the other part the search for all the compatible features or invariants in the χ 2 sense. We have observed that hashing techniques raise some important problems because of measurement uncertainty and of the non-Euclidean topology of the concerned manifolds. These problems turn crucial with invariants of features with high dimension. There is perhaps a trade-off to be found between adaptative space sampling and hashing methods.

In conclusion, generic, correct, and robust recognition algorithms on uncertain geometric features are possible. However, from a computational point of view, we need to design speed-up techniques for efficiency.

Figure 1 .

 1 Figure 1. Left:The tangent planes at points x and y of the sphere S 2 are different: the vectors v and w of T x M cannot be compared to the vectors t and u of T y M. Thus, it is natural to define the dot product on each tangent plane. Right: The geodesics starting at x are straight lines in the exponential map, and the distance along them is conserved.

VIDERE 1: 2

 2

VIDERE 1: 2

 2 Toward a Generic Framework for Recognition 63

VIDERE 1: 2

 2

VIDERE 1: 2

 2

VIDERE 1: 2

 2 Toward a Generic Framework for Recognition 68

Figure 2 .

 2 Figure2. Maximal and mean observed error for 3-D points due to the uncertainty of the trihedron in the basis for geometric hashing. Each of the 500 values is computed as the mean and maximal error for 1000 perturbations of a three-point basis. The x value is the corresponding predicted maximal error.

VIDERE 1: 2

 2

VIDERE 1: 2

 2 Toward a Generic Framework for Recognition 72

VIDERE 1: 2

 2 Toward a Generic Framework for Recognition 73

Figure 3 .

 3 Figure 3. Voronoï diagram and dual Delaunay triangulation on the sphere S 2 , generated interactively and in real time with the java applet ModeMap of Dave Watson (www.iinet.com.au/ ∼watson/modemap.html).

VIDERE 1: 2

 2 Toward a Generic Framework for Recognition 74

Figure 4 .

 4 Figure 4. Multidimensional data hashing with error.

Figure 5 .

 5 Figure 5. Left: a 512 • 512 (synthetic) 2-D image. The two points used for the basis are shown in bold with their error zone (five pixels). Right: the slice of the hash table corresponding to the measured length of the basis. One can see the sensitivity of the invariant coordinates by comparing the size of the error zones in the image space (left) and in the invariant coordinate space (right).

VIDERE 1: 2

 2

VIDERE 1: 2

 2

3

 3 Proportion of False Positives (Hough and Alignment)

VIDERE 1: 2

 2

VIDERE 1: 2

 2 Toward a Generic Framework for Recognition 81

Figure 7 .

 7 Figure 7. Qualitative estimation of the number of false positives involving at least x matches in MR images of 2,500 features. Comparison between frames and points: we need roughly five times more point matches than frame matches to obtain the same probability (10 frames and 56 point matches for a probability of 10 -10 ). The actual match found involves about 500 features, and its probability of being a false positive is thus practically zero.

VIDERE 1: 2

 2 Toward a Generic Framework for Recognition 82= DIAG(0.005, 0.006, 0.065; 0.45, 0.60, 0.165).

VIDERE 1: 2

 2 Toward a Generic Framework for Recognition 83

6 r • σ 6 t • π 3 24 = 4 . 8 ×

 62448 10 -2 for frames andV 0 = χ 3 • σ 3 t • 4•π 3 = 4.06. Finally, we obtain the following selectivities:η f r = 7.26 × 10 -7 and η pt = 1.2 × 10 -4 .

VIDERE 1: 2

 2

  Toward a Generic Framework for Recognition 62 and a trihedron (t 1 , t 2 , n), where (t 1 , t 2 ) ≡ (-t 1 , -t 2 ) are the principal directions of the surface.

The first two are not even rotations unless they are renormalized: the sum of orthogonal matrices is generally not an orthogonal matrix and the sum of unit quaternions is not a unit quaternion. Moreover, choosing the sign of the unit quaternion used in the sum is not a trivial problem. VIDERE 1:2

A frame is a point with an orthonormal trihedron and is equivalent to a Euclidean coordinate system.VIDERE 1:2Toward a Generic Framework for Recognition 67