Amaury Dame
email: eric.marchand@irisa.fr

Second order optimization of mutual information for real-time image registration Amaury Dame, Eric Marchand

Abstract-In this paper we present a direct image registration approach that uses Mutual Information (MI) as a metric for alignment. The proposed approach is robust, real-time and gives an accurate estimation of a set of 2D motion parameters. MI is a measure of the quantity of information shared by signals.

Although it has the ability to perform robust alignment with illumination changes, multi-modality and partial occlusions, few works propose MI-based applications related to spatio-temporal image registration or object tracking in image sequences due to some optimization problems that we will explain.

In this work, we propose a new optimization method that is adapted to the MI cost function and gives a practical solution for real time tracking. We show that by refining the computation of the Hessian matrix and using a specific optimization approach, the registration results are far more robust and accurate than the existing solutions while the computation is cheaper. A new approach is also proposed to speed up the computation of the derivatives and keep an equivalent optimization efficiency. To validate the advantages of the proposed approach, several experiments are performed.

Index Terms-Mutual information, registration, tracking, optimization.

I. INTRODUCTION

I Mage registration goal is to geometrically align two im- ages acquired at different time and from different camera viewpoints [START_REF] Brown | A survey of image registration techniques[END_REF], [START_REF] Zitová | Image registration methods: a survey[END_REF]. Considering a displacement model, this alignment process requires the optimization of a similarity measure. Various registration problems can be considered. First, one can consider the case where a wide baseline between two viewpoints is available. In this case, most of the approaches consist of the following steps: features or landmarks detection, features matching, displacement/transformation estimation. Possible applications for such registration methods include stereo-mapping to recover depth from disparities, remote sensing, mosaicing of a large area, medical image registration, etc. A good survey of such techniques is proposed in [START_REF] Zitová | Image registration methods: a survey[END_REF]. The second group of registration problem, also known as tracking, considers image sequence analysis where images differ only slightly and assumptions about smooth changes are justified. Although wide baseline registration techniques still apply, since a continuous motion is assumed from frameto-frame, other methods can be proposed and only a small increment of the transformation parameters has then to be estimated. Possible applications are motion estimation, video mosaicing, augmented reality, etc.

In this paper, we shall consider only registration in image sequences. Such approach which can be seen as a 2D motion estimation issue is also often referred as direct tracking or region tracking methods. Major difficulties in such a registration process are image noise, illumination changes and occlusions. Along with robustness to such perturbations, our motivation is to focus on registration and tracking considering different sensor modalities (eg, infra-red and visible images). The choice of a robust similarity measure is then fundamental. In this paper a process based on mutual information [START_REF] Shannon | A mathematical theory of communication[END_REF], [START_REF] Viola | Alignment by maximization of mutual information[END_REF] is proposed.

Most of the available direct tracking techniques can be divided into two main classes: feature-based and model-based registration method. The former approach focuses on tracking 2D features such as geometrical primitives (point, segments, circles, etc.) or object contours (such as active contours). The latter explicitly uses a model of the scene. This model can be a 3D model leading to a pose estimation process defined as a registration between the measures in the image and the forward projection of the 3D model [START_REF] Drummond | Real-time visual tracking of complex structures[END_REF] [START_REF] Comport | Real-time markerless tracking for augmented reality: the virtual visual servoing framework[END_REF]. One can also consider 2D models. Within this category, the features to be tracked can be represented by a descriptor. These descriptors can be image histograms leading to mean shift like approaches [START_REF] Comaniciu | Real-time tracking of non-rigid objects using mean shift[END_REF] or point neighborhood leading to keypoint tracking by matching approaches [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] [START_REF] Lepetit | Keypoint recognition using randomized trees[END_REF]. Following a statistical approach, [START_REF] Bibby | Robust real-time visual tracking using pixel-wise posteriors[END_REF] proposes an approach that merges both a level set approach and histogram-based approach to solve the registration problem. While very robust, these approaches are nevertheless not prone to the estimation of complex movements. In this attempt, it is possible to consider that the 2D model is a reference image (or a template). In that case, the goal is to estimate the motion (or warp) between the current image and a reference template. An example of such approaches are differential image registration methods such as the KLT [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] or its sequels [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] [START_REF] Hager | Efficient region tracking with parametric models of geometry and illumination[END_REF][1] [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF][20] [START_REF] Benhimane | Homography-based 2d visual tracking and servoing[END_REF]. Those approaches are not limited to 2D motion estimation, considering for example the motion of a planar object in the image, it is indeed possible to retrieve its 3D motion. The approach described in this paper is related to this last category of registration methods.

In the context of "KLT-like approaches", a measure of the alignment between the reference image and the current image and its derivatives with respect to the motion (warp) parameters is used within a non-linear estimation process to estimate the current object motion. What seems to be a well adapted measure is the standard Sum of Squared Differences (SSD) function [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] [START_REF] Baker | Equivalence and efficiency of image alignment algorithms[END_REF]. But such approaches are not effective in the case of illumination changes and occlusions. Several solutions have been proposed to add robustness toward those variations. Some include the use of M-estimators to deal with occlusions [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] [START_REF] Hager | Efficient region tracking with parametric models of geometry and illumination[END_REF] or add new parameters to estimate the illumination variations [START_REF] Hager | Efficient region tracking with parametric models of geometry and illumination[END_REF] [START_REF] Silveira | Real-time visual tracking under arbitrary illumination changes[END_REF]. One can also consider local normalized cross correlation (NCC) [START_REF] Irani | Robust multi-sensor image alignment[END_REF] or ZNCC to replace SSD.

In this paper, our goal is: first to have an image registration approach that is robust to environmental variations and second, that can handle multi-modality. The proposed solution is then no longer to minimize the SSD but to consider a more robust alignment function, the Mutual Information (MI) between the reference image and the current image, that is defined by the information shared by them and maximize it. MI has been introduced in the context of information theory by Shannon [START_REF] Shannon | A mathematical theory of communication[END_REF]. It has been later considered as an image similarity measure back in the mid ninety's independently by Collignon [9] for tomographic image registration, Studholme [START_REF] Studholme | Automated 3D registration of truncated mr and ct images of the head[END_REF] for MR and CT image, an by Viola [START_REF] Viola | Alignment by maximization of mutual information[END_REF] for projection image. Since then MI has become a classical similarity measure especially for multi-modal registration techniques [START_REF] Pluim | Mutual-information-based registration of medical images: a survey[END_REF] (eg, for medical or remote sensing applications). MI has proved to be robust to occlusions and illumination variations and, therefore, it can be also considered as a good alignment measure for image tracking. Nevertheless, to date only few works features mutual information within a frame-to-frame tracking method [START_REF] Chen | On registration of regions of interest (roi) in video sequences[END_REF], [START_REF] Dowson | A unifying framework for mutual information methods for use in non-linear optimisation[END_REF], [START_REF] Dowson | Mutual information for lucas-kanade tracking (milk): An inverse compositional formulation[END_REF], [START_REF] Panin | Mutual information-based 3d object tracking[END_REF]. To be used efficiently within a template tracking algorithm (such as [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF]), an optimization technique has to be considered. Various approaches have been proposed: first-order gradient descent [START_REF] Viola | Alignment by maximization of mutual information[END_REF], multi-resolution hill climbing algorithm [START_REF] Studholme | Automated 3D registration of truncated mr and ct images of the head[END_REF] or simulated annealing techniques [START_REF] Ritter | Registration of stereo and temporal images of the retina[END_REF]. Powell's [START_REF] Collignon | Automated multi-modality image registration based on information theory[END_REF], [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF] or Simplex [START_REF] Meyer | Demonstration of accuracy and clinical versatility of mutual information for automatic multimodality image fusion using affine and thin-plate spline warped geometric deformations[END_REF], [START_REF] Chen | On registration of regions of interest (roi) in video sequences[END_REF] methods (which do not require function derivatives to be analytically expressed) have been very popular in MI optimization but the former is sensitive to local optima in the registration criterion whereas the latter is known to be computationally inefficient. Considering that MI computation is evaluated from the joint image intensity histogram, an analytic derivative of the mutual information is difficult to obtain. In order to compute MI derivatives, [START_REF] Maes | Comparative evaluation of multiresolution optimization strategies for multimodality image registration by maximization of mutual information[END_REF] introduces partial volume interpolation for the construction of the joint histogram leading to an analytic computation of MI gradients. In [START_REF] Thévenaz | Optimization of Mutual Information for Multiresolution Image Registration[END_REF], the authors formulate the mutual-information criterion as a continuous and differentiable function of the registration parameters using B-Spline Parzen windows. These derivatives are then used in a Levenberg-Marquardt like minimization method. Such a formulation has been considered within a motion estimation process (or camera tracking) in [START_REF] Dowson | A unifying framework for mutual information methods for use in non-linear optimisation[END_REF], [START_REF] Dowson | Mutual information for lucas-kanade tracking (milk): An inverse compositional formulation[END_REF], [START_REF] Panin | Mutual information-based 3d object tracking[END_REF]. However the existing approaches are not taking full advantage of the accuracy of MI nor are they appropriate for real-time applications especially if a complex motion model is considered.

In this paper we present a MI-based image alignment. An important contribution is to propose an optimization process adapted to the MI cost function. We propose an inverse compositional optimization approach [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF] where an important part of the required derivatives can be precomputed, resulting in small computation times. A precise, complete and efficient computation of the Hessian matrix is described. We show that the inverse compositional approach allows the estimation of the Hessian matrix after convergence that can be used in a Newton's like approach to give an accurate and fast estimation of the displacement parameters. Finally a new approach is proposed to speed up the computation of the derivatives through a selection of the reference pixels making the image alignment process possible at video-rate.

In the remainder of this paper, Section II presents an overview of the differential image registration approaches. In section III, a brief introduction on information theory is given with the definition of mutual information, then a formulation adapted to the differential alignment method is presented. Section IV deals with the optimization of the resulting mutual information function with respect to the motion parameters to estimate. Finally section V presents several image registration or template tracking experiments including the Metaio benchmark and demonstrates the multi-modal capability of the approach.

II. DIFFERENTIAL TEMPLATE-BASED IMAGE

REGISTRATION

Differential image alignment [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF] (or template tracking) is a class of approaches based on the optimization of an image registration function.

The goal is to estimate the displacement p of an image template I * in a sequence of images. Considering a frame-toframe tracking process, the template I * is usually a region of interest extracted from the very first image of the sequence. In the case of a similarity function f , the problem can be written as :

p t = arg max p f (I * , w(I t , p)).
(

) 1
where we search the displacement p t that maximizes the similarity between the template I * and the warped current image I t . In the case of a dissimilarity function the problem would be simply inverted in the sense that we would search the minimum of the function f . For the purpose of clarity, the warping function w is here used in an abuse of notation to define the overall transformation of the image I by the parameters p. Afterwards, its proper formulation will be preferred using w(x, p) to denote the position of the point x transformed using the parameter p. The displacement parameters p can be of high dimension. For instance, the experiments that will be presented at the end of the paper consider a homography transformation that corresponds to p ∈ sl(3) that is 8 parameters. Approaches such as an exhaustive search of p are thus too expensive if not impossible.

To solve the maximization problem, the assumption made in the differential image registration approaches is that the displacement of the object between two consecutive frames is quite small. The previous estimated displacement p t-1 can therefore be used as first estimation of the current displacement to perform the optimization of f and incrementally reach the best estimation p t .

Multiple solutions exists to compute the update of the current displacement parameters and perform the optimization. Baker and Matthews showed that two formulations were equivalent [START_REF] Baker | Equivalence and efficiency of image alignment algorithms[END_REF] depending on whether the update is acting on the current image or the reference. The former is the direct compositional formulation which considers that the update is applied to the current image, thus we search the update ∆p that maximize f as:

∆p k = arg max ∆p f (I * , w(w(I t , ∆p), p k )). ( 2
)
This equation is typically solved using a Taylor expansion where the update is computed with the function derivatives with respect to ∆p. For a pixel x, the update of the current parameters p k is then applied as follows:

w( w(x, ∆p), p k ) → w(x, p k+1 ). (3)
A second equivalent formulation is the inverse compositional formulation which considers that the update modifies the reference image, so that ∆p is chosen to maximize:

∆p k = arg max ∆p f (w(I * , ∆p), w(I t , p k )). (4)
In this case the current parameters will be updated using:

w( w -1 (x, ∆p k ), p k ) → w(x, p k+1 ). (5)
In the inverse compositional formulation, since the update parameters are applied to the reference image, the derivatives with respect to the displacement parameters are computed using the gradient of the reference image. Thus, these derivatives can be partially precomputed and the algorithm is far less time consuming. Since we are interested in a fast estimation of the displacement parameters, the remainder of the paper will focus on the later inverse compositional approach. One essential choice remains the one of the alignment function f . One natural solution is to choose the function f as the sum of squared differences (SSD) of the pixel intensities between the reference image and the transformed current image:

p t = arg min p (SSD(I * , w(I t , p))) (6)
= arg min

p x∈W (I * (x) -I t (w(x, p))) 2 (7)
where the summation is computed on each point x of the reference template that is the region of interest (W) of the reference image. As suggested by its definition, this dissimilarity function is very sensitive to occlusions and illumination variations. Many solutions have been proposed to robustify the SSD. M-estimators can robustify the least squared problem toward occlusions [START_REF] Hager | Efficient region tracking with parametric models of geometry and illumination[END_REF], while a model of illumination changes can be coupled with the motion model to create a tracker robust to lighting changes [START_REF] Silveira | Real-time visual tracking under arbitrary illumination changes[END_REF]. Nevertheless those solutions are complex since additional parameters have to be estimated and aligning two images acquired using different modalities of acquisition remains impossible. Let us for example consider an aerial image and a map template (see figure 1(a)). Considering these two modalities is obviously an extreme case, but it will emphasize the robustness of the MI with respect to other similarity function. The value of SSD is computed with respect to the translations between the map and the satellite image. The two images are showing the same place (at least for a human eye they contain the same "information"), however, since the link between the intensities of the pixels is not linear, the SSD function represented in figure 1(b) gives no information on the alignment between the two images. Since NCC and ZNCC have shown some very good results in multi-modal alignment problems [START_REF] Irani | Robust multi-sensor image alignment[END_REF], we also evaluated the ZNCC efficiency in that matter. We can see in figure 1(c) that the case is too extreme and that ZNCC has also no significant optimum and therefore can not be used in this registration task.

To deal with occlusions, illumination variations and multimodality, we propose to define our alignment function as the mutual information [START_REF] Shannon | A mathematical theory of communication[END_REF], [START_REF] Viola | Alignment by maximization of mutual information[END_REF]. Originating from the information theory, MI is a measure of statistical dependency between two signals (or two images in our case) that is, as we will see, robust to all this variations of appearance.

III. MUTUAL INFORMATION

Rather than comparing intensities, mutual information is the quantity of information shared between two random variables. Mutual information of two random variable I and I * is then given by the following equation [START_REF] Shannon | A mathematical theory of communication[END_REF]:

M I(I, I * ) = h(I) + h(I * ) -h(I, I * ). (8)
where the entropy h(I) is a measure of variability of a random variable I (signal, image...). If r are the possible values of I and p I (r) = P (I = r) is the probability distribution function of r, then the Shannon entropy h(I) of a discrete variable I is given by the following expression:

h(I) = - r p I (r) log (p I (r)) . (9)
The probability distribution function of the gray-level values is then simply given by the normalized histogram of the image I. The entropy can therefore be considered as a dispersion measure of the image histogram. Following the same principle, the joint entropy h(I, I * ) of two random variables I and I * can be defined as the variability of the couple of variables (I, I * ). The Shannon joint entropy expression is given by: notation as in section II, the mutual information can thus be written with respect to p:

h(I, I * ) = - r,t p II * (r, t) log (p II * (r, t)) (10
M I(p) = M I(w(I, p), I * ) = h(w(I, p)) + h(I * ) -h(w(I, p), I * ). ( 11
)
The final expression of MI is obtained by developing the previous equation using the entropy equations ( 9) and ( 10):

M I(p) = r,t p II * (r, t, p) log p II * (r, t, p) p I (r, p)p I * (t) (12)
The analytical formulation of a normalized histogram of an image I * is classically written as follows:

p I * (t) = 1 N x x φ t -I * (x) (13)
p I (r, p) = 1 N x x φ r -I(w(x, p)) (14)
p II * (r, t, p) = 1 N x x φ r -I(w(x, p)) φ t -I * (x)
where x are the points of the region of interest in the image and N x is the number of points. r and t are the possible values of I(x) and I * (x) that are the scaled version of the original images, so that {r, t} ∈ [0, N c ] 2 . Let us note that to have a smooth mutual information it is important to maintain the number of histogram bins N c low (N c = 8 in our implementation). In the classical formulation φ is a Kronecker's function: φ(x) = 1 for x = 0 and φ(x) = 0 otherwise. So that, for instance, each time I * (x) = i the ith histogram bin value of p I * is incremented. However this formulation does not take advantage of the decimal part of the scaled intensities, therefore several solutions have been proposed to simultaneously smooth the mutual information function, make its formulation differentiable and keep its accuracy [39][25]. Several approaches propose to use Gaussian function, however, in our approach, we focus on the use of B-splines functions for φ [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF], [START_REF] Thévenaz | Optimization of Mutual Information for Multiresolution Image Registration[END_REF]. It has indeed been shown that these functions provides a good approximation of Gaussian functions while their computation and the one of their derivatives is cheaper. As it will be discussed later, this also permits to have a smooth, accurate but computationally cheap gradient of MI to perform its optimization.

IV. MUTUAL INFORMATION-BASED MOTION ESTIMATION

In this section we will see how to use the MI cost function with the differential image registration formulation presented in section II. Once our approach is fully defined, a pseudo-code of the algorithm is given to summarize the proposed method.

A. Overview

The goal of our tracking problem is to align an image template I * with an input image I. If we assume that the reference template appears in I, the goal is to search for the transformation that aligns the pixels x of the reference image I * to the corresponding pixels x of I in the sense of our chosen similarity measure. Assuming that the transformation from the reference points to the input image can be modeled by a warp function x = w(x, p), the problem can be formulated as:

p = arg max p M I (I * (x), I(w(x, p))) (15)
Since this problem is impossible to solve linearly, a non-linear optimization is performed. To initialize the optimization, a first guess of the displacement parameters is required. Since we suppose that the displacement of the object between two consecutive frames is small, a good approximation is to approximate the parameter p t of the input image I t at a time t using the parameters estimated for the previous frame:

p t = p t-1 .
To initialize the whole tracking approach, the position of the template in the first image I 0 has to be known coarsely. Since the first image of the sequence is usually the one that defines the template I * , the first displacement parameters p 0 between the template and the first image simply correspond to an identity transformation (considering the warp functions defined in the first chapter it yields: p 0 = 0). Otherwise, the first estimation can be performed using some matching process, such as a keypoints matching approach [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], [START_REF] Lepetit | Keypoint recognition using randomized trees[END_REF] or other wide baseline registration method. The first approximation of the displacement p 0 t = p t-1 is then refined using the numerical resolution of the equation [START_REF] Drummond | Real-time visual tracking of complex structures[END_REF]. To solve the maximization, an iterative optimization method is used that successively goes closer and closer to the optimum of the cost function p t (see Figure 2). For a clarity purpose, let us now consider the maximization peculiar to one image I (we drop the subscript t) and focus on the iteration number noted using the superscript k.

Let us recall that for efficiency issue, we chose to consider an inverse compositional approach. The difference with the forward compositional approach comes from the optimization process where the updating steps from the current guess to

Previous image I t-1
Current image I t Fig. 2. The first approximation of the position p 0 t , given by the previous position p t-1 , is iteratively refined to find the optimal parameters pt.

the optimal displacement parameters are modified. Instead of searching the update parameters that will bring the warped points of the current image into the points of the template image, the formulation of the problem is inverted so that we search the "inverse" update that brings the points of the template image into the warped points of the current image. In the inverse compositional approach [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF], the goal is then formulated as finding the update ∆p that leads to the optimum, so that, at each iteration k, we seek (see Section II for details):

∆p k = arg max ∆p M I I * (w(x, ∆p)), I(w(x, p k )) (16)
The optimization using this formulation is similar to the optimization using the forward compositional approach. Nevertheless, since the update is considered to affect the reference image, we will see that more elements of the mutual information derivatives with respect to the update can be precomputed.

B. Derivative function analysis

Let us remind that the goal is to estimate the displacement parameters p t that maximizes the MI using a first estimation of the parameters p t-1 and an iterative update of the parameters. In this work we ought to register planar regions through 3D displacements. This problem implies a strong correlation between the elements of the vector p. Therefore, estimating the update using first-order optimization method such as a steepest gradient descent is not adapted. Such non-linear optimization are usually performed using a Newton's method that assume the shape of the function to be parabolic.

Newton's method uses a second order Taylor expansion at the current position p k-1 to estimate the update ∆p required to reach the optimum of the function (where the gradient of the function is null). The same estimation and update are performed until the parameter p k effectively reaches the optimum. The update is estimated following the equation:

∆p = -H -1 G ( 17
)
where G and H are respectively the gradient and Hessian matrices of the mutual information with respect to the update ∆p. Following the inverse compositional formulation defined in equation ( 4) those matrices are equal to:

G = ∂M I(w(I * , ∆p), w(I, p)) ∂∆p (18)
H = ∂ 2 M I(w(I * , ∆p), w(I, p)) ∂∆p 2 (19)
Applying the derivative chain rules to equation ( 12) yields the following gradient and Hessian matrices:

G = r,t ∂p II * ∂∆p 1 + log p II * p I * (20) H = r,t ∂p II * ∂∆p ∂p II * ∂∆p 1 p II * - 1 p I * + ∂ 2 p II * ∂∆p 2 1 + log p II * p I * (21)
For the purpose of clarity, the marginal probabilities and joint probability that are actually depending on r, t, p * and ∆p are simply denoted as p I , p I * and p II * . The details of the calculation from equation [START_REF] Irani | Robust multi-sensor image alignment[END_REF] to equation ( 21) can be found in [START_REF] Dowson | A unifying framework for mutual information methods for use in non-linear optimisation[END_REF].

By analogy with the Hessian computation in a Gauss-Newton's method for a least squared problem that is assuming that the neglected term is null after convergence, second order derivatives are usually neglected in the Hessian matrix computation [START_REF] Thévenaz | Optimization of Mutual Information for Multiresolution Image Registration[END_REF], [START_REF] Thévenaz | Optimization of mutual information for multiresolution image registration[END_REF], [START_REF] Dowson | A unifying framework for mutual information methods for use in non-linear optimisation[END_REF], [START_REF] Dowson | Mutual information for lucas-kanade tracking (milk): An inverse compositional formulation[END_REF] leading to:

H r,t ∂p II * ∂∆p ∂p II * ∂∆p 1 p II * - 1 p I * . ( 22
)
In our approach we compute the Hessian matrix using the second order derivatives. In our point of view, they are required to obtain a precise estimation of the motion. Indeed, let us consider the approximation made in [START_REF] Lieberknecht | A dataset and evaluation methodology for template-based tracking algorithms[END_REF]. Considering the expression of the marginal probability p I * (t) = r p II * (r, t), it is clear that p I * (t) > p II * (r, t) so 1/p II * (r, t)-1/p I * (t) > 0.

Since ∂p II * ∂∆p ∂p II *
∂∆p is a positive matrix then the final Hessian matrix given by ( 22) is positive. Since the goal is to maximize MI, the Hessian matrix after convergence is supposed to be negative by definition. The common approximation of ( 22) is thus not suited for the optimization of MI.

As we can see in equation [START_REF] Jurie | Hyperplane approximation for template matching[END_REF] and equation ( 21), the derivatives of the mutual information depend on the derivatives of the joint probability. Using the previous definition in [START_REF] Dowson | Mutual information for lucas-kanade tracking (milk): An inverse compositional formulation[END_REF] and passing the derivative operator through the summation yields the following expressions:

∂pII * ∂∆p = 1 Nx x φ t -I(w(x, p)) ∂φ r -I * (w(x, ∆p)) ∂∆p ∂ 2 pII * ∂∆p 2 = 1 Nx x φ t -I(w(x, p)) ∂ 2 φ r -I * (w(x, ∆p)) ∂∆p 2 . ( 23
)
The remaining expressions to evaluate are the variations of the B-spline function φ with respect to the update. Their derivatives are obtained using the chain rule leading to:

∂φ(r -I * (w(x, ∆p))) ∂∆p = - ∂φ ∂r ∂I * ∂∆p (24)
∂ 2 φ r -I * (w(x, ∆p))

∂∆p 2 = ∂ 2 φ ∂r 2
∂I * ∂∆p

∂I * ∂∆p - ∂φ ∂r ∂ 2 I * ∂∆p 2 . ( 25
)
Finally the derivatives of the reference image intensity with respect to the update parameters ∆p is given by the following expressions:

∂I * ∂∆p = ∇I * ∂w(x, p) ∂∆p (26)
∂ 2 I * ∂∆p 2 = ∂w ∂∆p ∇ 2 I * ∂w ∂∆p + ∇I * x ∂ 2 w x ∂∆p 2 + ∇I * y ∂ 2 w y ∂∆p 2 (27)
where ∇I * are the image gradients of the reference image, obtained using the convolution of a Gaussian filter and a derivative filter, the Gaussian filter allowing for a smoother version of the gradients. The motivation for using the inverse compositional formulation is then obvious. The derivatives of the warp function are all computed at ∆p = 0, their values are then constant for each pixels of the template. Moreover, since the reference image I * is constant, its gradients and all the expressions from equation [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] to equation ( 27) are constants and have to be precomputed only one time.

In our work we focus on planar region registration. The warp function is thus defined by the group action w : SL(3) × P 2 with x ∈ P 2 and p defines the 8 parameters of the sl(3) lie algebra associated to the SL(3) group. However, this research is not limited to such a warp function thus details will not be given on the warp derivatives. All details regarding the derivatives of the chosen warp function can be found in [START_REF] Benhimane | Homography-based 2d visual tracking and servoing[END_REF].

Let us emphasize that any kind of warp model can be considered. Although homography have been considered in this paper, it can also be applied on affine motion model [START_REF] Dame | A unified direct approach for visual servoing and visual tracking using mutual information[END_REF], pose parameters SE(3) [START_REF] Panin | Mutual information-based 3d object tracking[END_REF] [START_REF] Caron | L'information mutuelle pour l'estimation visuelle directe de pose[END_REF] and other motion models. The method could also be extended to non-rigid registration process. In that case, specific local distortions have to be considered. Radial basis functions (such as Wendland's function or thin-plate splines) are able to handle locally varying geometric distortions and can be considered within the proposed framework. In any cases, although the main algorithm will remain unchanged, the only modification will be to redefine the warp derivatives.

C. Optimization approach

The Newton's method that can be used to perform the estimation of the update parameters ∆p is based on the assumption of a similarity function with a parabolic shape. One can immediately notice that this assumption can be easily violated by looking at the function's shape where we see that the assumption is correct only near the maximum. Since the violation could cause the Newton's method to fail, a better approach has to be found.

To evaluate the efficiency of the following optimization methods, a set of alignment experiments has been realized. The goal is to estimate the known position p * of a template in an image (see figure 4(a)) from many initial position parameters (see figure 4(b)). The initial parameters are automatically generated applying a random noise to the ground truth position.

The convergence rate of the optimization method are then evaluated with respect to the initial positioning error. The positioning error err is defined as the RMS distance between the correct position of some reference points

x * i = w(x i , p * )
and the current position of the points w(x i , p) [START_REF] Lieberknecht | A dataset and evaluation methodology for template-based tracking algorithms[END_REF]. The reference points are simply chosen as the 4 corners of the template so that the error becomes:

err(p) = 4 i=1 x * i -w(x i , p) (28)
We consider that the optimization converges as soon as the error err is below 0.5 px. 500 alignment experiments are performed for each initial positioning error err from 1 to 20 that is a total of 10000 experiments. As output we retrieve the convergence rate, the average number of iterations required to reach convergence, the final residues after convergence and the computation time of each iteration. Indeed, those values gives a good overview of the efficiency of the optimization methods.

The Gradient descent method cannot estimate an accurate estimation of the homography (see section IV-B). Indeed its use gives a final estimation with an error always above 0.5 px for the all set of experiments (that is a 0% convergence rate). Thus the results have not been included in figure 4.

1) Newton's method: Mutual information function is a quasi-concave function, thus the parabolic hypothesis of the Newton's method is only valid near the convergence. As soon as the displacement in the sequence is important, the initial parameters p t-1 would be on the convex part of the cost function that will cause the optimization to diverge.

The problem is in fact equivalent using a SSD function. One example of the values obtained on the estimation of a translational displacement is presented in figure 3 for both the MI function and the minus of the SSD function. For the purpose of clarity, we choose to analyze the minus of the SSD function to deal with a maximization for both functions. The quasi-concave shape of both functions is obvious. The parabolic assumption is only correct for the concave part of the function, that is where their second order derivatives are negative (the area highlighted in purple). Therefore the convergence domain using a classical Newton's method would be very small. Figure 4(c) shows the convergence results obtained using our set of convergence tests. The convergence domain of the Newton's method is indeed practically very small in the case of the homography estimation. As soon as the initial error exceeds 2 px, the initial parameters are, most of the time, out of the convergence domain of the Newton's method and the convergence rate decreases drastically. Considering the converging experiments, then once the convergence is achieved, the parabolic shape assumption is verified and the method gives good quality estimation with a mean final residue of 0.06 px. However, it is rarely the case and the computation of the Hessian at each iteration makes the process really time consuming (see figure 5 for the time per iteration).

Considering the simple one dimensional example, one could expect an optimization that has a convergence domain as wide as the one of the gradient descent method (the blue area in figure 3).

2) Conditioning the optimization: In this section, we show how to combine the convergence domain of the gradient Fig. 3. SSD, MI and their derivatives with respect to one translation (px). The purple area is the convergence domain using a classical Newton's method, the blue one is the convergence domain of a Gradient descent method. The proposed method keeps the wider convergence domain of the gradient's method in blue, while having the convergence properties of the Newton's method near the optimum, allowing an accurate estimation of complex transformations.

Ground truth

Initial positions (err = 20px) descent with the accuracy and efficiency of the Newton's method for the optimization of MI. In registration problems formulated with a SSD function, the Gauss-Newton approximation condition the problem by estimating a Hessian matrix that is always definite positive (see section IV-B and the green curve in figure 3) and that is a good approximation of the exact Hessian matrix after convergence. Therefore its use permits to have a convergence domain as wide as the one with a gradient method (blue area) and a good convergence behavior next to the optimum yielding to accurate estimations.

In the mutual information maximization, the problem is different. Indeed approximating the Hessian matrix as it is proposed in [START_REF] Thévenaz | Optimization of Mutual Information for Multiresolution Image Registration[END_REF], [START_REF] Dowson | A unifying framework for mutual information methods for use in non-linear optimisation[END_REF], [START_REF] Dowson | Mutual information for lucas-kanade tracking (milk): An inverse compositional formulation[END_REF] do not gives an estimation of the Hessian matrix after convergence (see the green curve in 3 for the MI function). No approximation on the Hessian of MI simplifies the problem as the Gauss-Newton approach does for the SSD.

The solution that we propose is inspired from the Gauss-Newton approach. The idea remains to use an estimation of the Hessian matrix after convergence. To compute this estimation we consider that after convergence the alignment between the template and the warped current image is perfect. Therefore we simply assume that at the optimum we have I(w(x, p)) = I * (x).

This solution has several advantages:

• It gives a definite negative Hessian matrix that yields to have a wide convergence domain (blue area in figure 3). We can notice that the resulting convergence domain is as wide as the one of the SSD function in the considered 1D example. In section V-A2, further experiments will show that it is also the case for a homography estimation. • Since the Hessian matrix used in the Newton's method is the Hessian matrix after convergence, the behavior of the optimization near convergence is optimal and the final estimated displacement parameters are very accurate. • This approach has the advantage of its computation time.

In the classical Newton's method the Hessian and Jacobian are computed for each iterations. In the proposed approach the Hessian matrix is computed one time in the whole experiment. The proposed optimization has been evaluated on the set of experiment presented in figure 4. As expected, the convergence domain is larger than the one using the classical Newton's method. The optimization converges for all the experiments with an initial error below 16px and the convergence rate slightly decreases for err > 16. Since we use the Hessian estimated after convergence, then the behavior near the optimum is suited to reach an accurate solution yielding to final residues that have a mean of 0.06 px.

Figure 4(d) shows the number of iterations to reach convergence. The number of iterations with the proposed method is fewer than the one with the classical Newton's method while its computation is much cheaper (see figure 5).

3) Improving the computation time: Compared to a simple least squared problem, mutual information can still be considered as a very complex function to compute. The proposed approach offers already a practical solution. Nevertheless, faster performance is sometimes desired.

To compute the MI between the two images, all the information is required, so all the reference pixels must be used to compute the marginal and joint probabilities. As for the variation of the mutual information computation, only the motion of the pixels that are not in a uniform region will have a strong effect. This fact is obvious from equation ( 26) and [START_REF] Meyer | Demonstration of accuracy and clinical versatility of mutual information for automatic multimodality image fusion using affine and thin-plate spline warped geometric deformations[END_REF]. One very simple modification is then to perform the computation of the gradient and Hessian using only a selection of pixels in the template.

A simple measure to determine if a point is in a uniform region of the template is given by the norm of the reference image gradients. Therefore the selection condition can be written as:

∇I * (x) > α ( 29
)
where α is a given threshold. The summation in equation ( 23) is therefore computed on the reference pixels that respect this condition.

The efficiency of the proposed approach has been compared to the previous one using the set of experiments represented in figure 4. Using a threshold α = 6, the selected number of points corresponds to 18% of the total number of reference points. We can see on figure 4(c &d) that the convergence rate and the number of required iterations is equal to the ones of the previous method up to few percent and iterations. Curves "proposed" and "fast proposed" are superimposed meaning precision efficiency remains the same (the final residues still has a mean value of 0.06 px) whereas computational efficiency greatly improves: 3.5ms per iteration vs 5.1ms (see figure 5) for the same number of iterations (see figure 4(d) ).

In summary, for a similar efficiency, the computation time of the proposed method is 30% smaller (see figure 5). Such a selection method is therefore highly recommended in MI derivatives computation.

To summarize the whole process, a pseudocode of the algorithm is presented in Figure 6 in the case of a classical image registration task.

V. EXPERIMENTAL RESULTS

The differential image registration method that is presented in this paper has been implemented on a laptop with a 2.4GHz processor. The evaluation of the displacement parameters has been performed using the presented inverse compositional scheme combined with a pyramidal approach that increases the convergence domain and speeds up the convergence of the optimization. We limit our experiments to the estimation of the displacement of planar image regions.

A. Mono-modal image alignment

The robustness and accuracy of the proposed method have been evaluated on various image sequences.

1) Image alignment through natural variations: This experiment concerns an indoor sequence acquired at video rate (25Hz). The initialization of the registration process has been performed by learning the reference image from the first image of the sequence and setting the initial homography to an identity. The template includes 16000 reference pixels.

The sequence has been chosen to illustrate the robustness of the motion estimation through many perturbation. Some images of the sequence are shown in figure 9. Firstly, the object is subject to several illumination variations: the artificial light produced an oscillation on the global illumination of the captured sequence. Moreover the object is not Lambertian, thus the sequence is subject to saturation and specularities (see figure 9 frame 200). The object is moved from its initial position using wide angle and wide range motions (figure 9 frame 400). And finally the object is subject to fast motion causing a significant blur in many images (figure 9 frame 600).

The frames of the sequence are presented with the corresponding estimated positions of the reference image. No ground truth of the object position is known, however, the projection of the tracked image on the reference image has been performed and qualitatively attests the accuracy of the registration process. Indeed the reconstructed templates show strong variations in terms of appearance but not in terms of position. We can conclude that the estimation of the motion is robust and accurate despite the strong illumination variations and blurring effects.

Concerning the processing time, using the proposed approach with no selection of the reference points (section IV-C2), the images are processed at video rate (25Hz). Using the fast computation (section IV-C3) it is about 40Hz. All the corresponding sequences are presented in the attached video.

In this experiment, we see that even if nothing guaranties that the optimization reaches the global maximum, the proposed computation of MI has such a wide maximum that it yields to a really robust approach. If nonetheless the convergence had to be verified, a solution could be to use a parallel tracking by matching approach [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] and check if we can find a better match (an higher MI score) than the one estimated in the non-linear optimization.

2) Evaluation on benchmark datasets: To have a quantitative measure of its accuracy and robustness, the registration process has been evaluated on some very demanding reference datasets proposed by Metaio GmbH [START_REF] Lieberknecht | A dataset and evaluation methodology for template-based tracking algorithms[END_REF]. Those datasets include a large set of sequences with the typical motions that we are suppose to face in augmented reality applications. Indeed sequences using eight reference images from low repetitive texture to highly repetitive texture are included. And for each reference image is a set of four sequences depicting wide angle, high range, fast far and fast close motion and one sequence with illumination variations.

The estimated motion has been compared with the ground truth for each sequences. The percentages given in the tables have been computed by Metaio relative to their ground truth. The upper table on figure 8 shows the results that have been obtained using the proposed approach. The tracker is considered as converging if the error between the estimation and the ground truth is below a given threshold. The error measure is similar to the one defined in equation [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]. A detailed definition is available in [START_REF] Lieberknecht | A dataset and evaluation methodology for template-based tracking algorithms[END_REF]. Some images of the sequences are shown in Figure 7 with the estimated position of the reference template. The mutual information based tracker proves its robustness and accuracy on most of the sequences.

The results obtained using the ESM approach [START_REF] Benhimane | Homography-based 2d visual tracking and servoing[END_REF] reported from [START_REF] Lieberknecht | A dataset and evaluation methodology for template-based tracking algorithms[END_REF] are also represented in the lower table of figure 8 where better convergence results are in bold characters. If we compare the results of the two methods we can see that both have similar convergence rates in most cases. But MI has an undeniable advantage in the cases of illumination variations

experiments.

From this demanding set of experiments, we can conclude that the proposed MI computation has a large convergence domain (at least as large as the one in the least squared problem) and that the proposed optimization is adapted to use the potential of the MI function leading to a very efficient image registration method.

B. Multi-modal image alignment 1) Satellite images versus map: This experiment illustrates the capabilities of the presented mutual information-based image registration process in alignment applications between map and aerial images. The reference image is a map template provided by IGN (Institut Géographique National) that can easily be linked to Geographic Information System (GIS) and the sequence has been acquired using a moving USB Fig. 8. Ratio of successfully registered images for our approach compared to the ESM [START_REF] Lieberknecht | A dataset and evaluation methodology for template-based tracking algorithms[END_REF].

camera focusing on a poster representing the satellite image corresponding to the map.

As it has been previously noticed in figure 1, a non-linear relationship exists between the intensities of the map and aerial image and this link can be evaluated by the MI function. Mutual information can therefore allow for aligning the satellite image using the map image. As Figure 10 shows, the selected initial position can be rather far from the correct position. Figure 11 shows the reference image and some images of the sequence with the corresponding overlaid results. There is no available ground truth for this experiment, nevertheless the overlaid results give a good overview of the alignment accuracy. We can also see in the attached video that the alignment converges despite some strong blurring effects.

2) Airborne infrared image versus satellite images: The same method has been evaluated with another current modality. This time the reference is a satellite image and the sequence is an airborne infrared sequence provided by Thales Optronic.

As expected, although very different, the two images shown in figure 13 are sharing a lot of information and thus MI can handle the registration process between the infrared sequence and the satellite image template. The warp function is still a homography. The airport scene is then supposed to be planar leading to an approximation. Nevertheless the proposed method remains robust. No ground truth is available, but the overlaid images qualitatively validates the accuracy of the registration process. The homographies have been decomposed to estimate the position of the plane with respect to the airport. The resulting 3D trajectory of the camera is represented in figure 12, as we can see the trajectory is smooth and has the expected behavior that shows the approach of a plane with respect to the runway. The trajectory of the camera with respect to the time is presented in the attached video. Figures 13 also shows some registered images that validate the accuracy of the motion estimation. The complete sequences are visible in the attached video.

To illustrate the improvements led by our approach Figure 14 shows the difference between a classical first order maximization approach using MI [START_REF] Dowson | A unifying framework for mutual information methods for use in non-linear optimisation[END_REF] (first row) and the proposed one (second row) that considers the full computation of the Hessian. Small registration errors can be observed when considering the classical approach while using a complete Hessian allows a better estimation of the transformation. Plots on the right of Figure 14 show the estimated altitude of the camera/plane (up to a scale factor) during the landing step. One can see, that in the first case, the estimation of the trajectory is noisy while with our approach one can clearly identify the classical three different steps of a landing process: airplane reduces downward slope, classically from 5 to 3 degree, this Reference template Alignment result Fig. 13. Registration of a satellite template image using MI on an airborne infrared sequence. 6 frames are represented with the over-imposed aerial reference (inside the green rectangle) projected using the estimated homography (Infrared images courtesy of Thales, optical image is obtained from google earth).

can be seen at 460 on the abscissa axis and finally executes a flare (350). Figure 13 shows the reference template (left), and six images of the sequence (overlaid with the projected template).

When estimating an homography (that is 8 independent parameters), template size may be important especially when considering multi-modal images. To sufficiently populate the joint probability histogram necessary for the computation of the mutual information and to maintain a good accuracy, we have always considered at least 100x100 pixels template. Nevertheless, although the computational cost will obviously increase with the number of pixels, it has to be remembered that most of the computation (Gradient, Hessian) are precomputed and the approach remains cost effective.

C. Mosaicing application

Image mosaics are a collection of overlapping images. The goal of the mosaicing problem is to find the transformations that relate the different image coordinates. Once the transformation between all the images is known, an image of the whole scene can be constructed. This problem requires to find a warping function that maps the coordinates of one image into the coordinate system of another image. When considering a video, one has to warp each new image into the coordinate system of the very first image of the video [START_REF] Szeliski | Image alignment and stitching: a tutorial[END_REF][6] [19][16]. This is basically a motion estimation process performed on the whole image. One can consider to estimate this motion using matched keypoints as in [START_REF] Brown | Automatic panoramic image stitching using invariant features[END_REF] or using SSD based motion estimation as in [START_REF] Irani | Mosaic based representations of video sequences and their applications[END_REF] [START_REF] Gracias | Underwater video mosaics as visual navigation maps[END_REF]. The latter approach is very efficient when image sequences are considered, that is, when displacements between one frame to another are small but shows its limits in case of image noise or occlusions. This section shows the benefits of using the presented approach to solve the registration problem.

These experiments show the application of the MI -based motion estimation algorithm to the mosaicing problem. In these sequences, since some parts of the scene completely disappear, it is necessary to define multiple reference images along the sequence. The approach is build as follows:

• Initialization: the first image is chosen as reference image, i.e. I * 0 = I 0 . • Registration: for every frame, we compute the displacement p k between I t and I * k . • Reference Update: every 30 images, the reference image I * k is changed and defined as the current image, i.e. I * k = I t for t = 30k (involving a small drift).

Using the homography from the current image to the current reference image and the homographies between the references, we retrieve the homography between the current image and the first image. Using this homography, we can project all the images of the sequence into the mosaic image and construct the global image of the whole scene.

In the first experiment (Figure 15 and 16), the overlapping images are simply a compressed sequence of 1000 images. The aerial scene is acquired from a camera embedded on a flying UAV and shows the ground that is approximately 500 meters away from the camera. Since this distance is very large, the scene can be approximated as a plane and registered using homographies. During the acquisition of the sequence, the camera is moving forward and is rotating around the vertical axis.

In Figure 15 we show some images from the sequence. This sequence has been downloaded on Youtube and is affected by the H264 coding artifacts. We can also note the poor quality of the images. Despite this poor quality, the resulting mosaic presented in Figure 16 shows the accuracy of the MI based method. Since the camera is making an entire revolution, the first and last images are overlapping. A very small drift occurs between the first and last estimated positions. Let us note that nothing has been performed to reduce the drift (such as the bundle adjustment approach proposed by [START_REF] Brown | Automatic panoramic image stitching using invariant features[END_REF]). Considering the template update problem and the planar assumption, the estimated homographies are accurate. A second similar experiment presents a mosaic build using more than 10000 images. The images are extracted from a highly compressed video. The camera was attached to a free flying balloon flying over Paris. Figure 17 shows three steps of the mosaic construction. In the last experiment (see mosaic in Figure 19), we consider Registration with a first order minimization approach using MI Registration with full computation of the Hessian Fig. 14. Registration of a satellite template image using MI on an airborne infrared sequence. First row corresponds to a registration process with a first order minimization approach while second row depicts our approach with a full computation of the Hessian. On the left is the estimated altitude versus distance to the runway (up to a scale factor). Fig. 19. Mosaic created from the John Ford movie "she wore a yellow ribbon". An affine motion model was considered. Note that some cavalrymen are moving all along the sequence. Despite these disturbances, motion is correctly estimated. a sequence extracted from the John Ford movie "she wore a yellow ribbon". To build this mosaic, an affine motion model was considered. The interest of this sequence is that some cavalrymen are moving all along the sequence and, therefore, act as important occlusions as can be seen on Figure 18. A comparison with a mosaic built using the SSD criterion is proposed in Figure 20 and demonstrates the robustness of our approach.

VI. CONCLUSION

This paper presents a robust and accurate template basedmotion estimation process that was defined using a new approach based on the mutual information alignment function. The definition of MI has been adapted to the differential image alignment problem so that the function is smooth and as concave as possible. The proposed definition preserves the advantages of MI with respect to its robustness toward occlusions, illumination variations and images from different modalities. A new optimization approach has been defined to deal with the quasi-concave shape of MI. The proposed approach is taking advantage of both the wide convergence domain of MI and its accurate maximum and besides is not computationally expensive. Moreover the time consumption is greatly reduced using a new approach based on the reference pixels selection that yields to an accurate, fast and robust registration process.

Finally the proposed method has been evaluated using several experiments. Its robustness and accuracy is verified using reference datasets and shows its advantages compared with classical approaches on monomodal image registration method. Some new applications are also proposed to use a model image acquired from another modality than the original sequence.

The algorithm presented here has been limited to planar scene. Nevertheless the proposed approach could similarly be applied to more complex model-based tracking applications where we could directly estimate the position of the object on SE(3) [7] [START_REF] Panin | Mutual information-based 3d object tracking[END_REF]. The method could also be extended to non-rigid registration process.

Fig. 1 .

 1 Fig. 1. Alignment functions wrt. translations between two images from the same area: (a) aerial image and the map reference. MI shows a maximum near zero translation at the alignment position whereas SSD and ZNCC gives no clear information on the alignment quality.

Fig. 4 .

 4 Fig. 4. Empirical convergence analysis of the optimization methods. The proposed methods (blue and green curves) have a very high convergence rate compared to the classical Newton's methods (red curve).

Fig. 5 .

 5 Fig. 5. Average time in ms required to compute one update ∆p of the optimization methods.

  %Definition of the reference: I * = I 0 ; p 0 = Id; %Precomputation of the derivatives: for x\{ ∇I * (x) > α} do Compute ∂φ(r-I * (w(x,∆p))) ∂∆p ; end for Compute H -1 using I(w(x, p)) = I * (x); for each image I t do p 0 t = p t-1 ; k = 0; while ∆p is significant do Compute G; ∆p = -H -1 G ; p k+1 t = p k t ⊕ ∆p -1 ; %See equation (5) end while end for

Fig. 6 .Fig. 7 .

 67 Fig. 6. Pseudo-code of the proposed method to solve a classical image registration task.

Fig. 9 .

 9 Fig. 9. registration of a planar template through illumination variations. First row: frame 0, 200, 400 and 600. The green rectangle represents the rectangle from the template image transformed using the estimated homography. Second row: projection of the templates for the same iterations in the reference image.

Fig. 10 .

 10 Fig. 10. Registration of an aerial sequence with a map template image by MI: (a) the considered template to be registered with an image of an aerial view of the same area (b) (image and map source: IGN). (c) Shows the initialization process (corresponding to one step of the tracking algorithm. It allows to show that the convergence domain is quite large despite the fact that the images are very different.

Fig. 11 .

 11 Fig. 11. Registration of an aerial sequence with a map template image by MI: frames 1, 250, 500 and 750 are represented with the over-imposed satellite reference (inside the green rectangle) projected using the estimated homography (image and map source: IGN) .

Fig. 12 .

 12 Fig.[START_REF] Dame | A unified direct approach for visual servoing and visual tracking using mutual information[END_REF]. From the homography to the estimation of the camera position. Green curve: estimated camera trajectory in the 3D space, blue: the 6 estimated camera positions corresponding to the frames represented in figure13.

Fig. 15 .

 15 Fig. 15. Some overlapping key images used for the mosaicing application.

Fig. 16 .

 16 Fig.[START_REF] Gracias | Underwater video mosaics as visual navigation maps[END_REF]. Resulting mosaic image: despite the poor quality of the sequence and the approximation that the scene is planar, the final displacement between the first and last image is accurate.

Fig. 17 .

 17 Fig. 17. Three steps of the "Paris" mosaic construction. The sequence feature more than 10000 images acquired from a camera attached to a free-flying balloon.

Fig. 18 .

 18 Fig. 18. Three images used for the "yellow ribbon" mosaic construction.

Fig. 20 .

 20 Fig. 20. Mosaic created from the John Ford movie but the similarity function is the SSD.

ACKNOWLEDGMENT

This work was supported by DGA under contribution to student grant.