*2.2. Automatic Registration*

The spatial mapping *T* used for the image morphing is determined by the image registration. Several registration methods are available. However, many of them require to define manually a set of corresponding points from the images *u* and *v*. We are interested in an automatic registration procedure that only needs the images *u* and *v* as inputs without any extra specifications. This requires the images to be similar enough for the automatic registration procedure to work.

We use the method described by Reference [33] based on the minimization of a cost function *J* with respect to the mapping *T*. The cost function can be divided in two terms (Equation (9)). The first one (*Jo*) represents the mapping error between the displaced original signal *u* ◦ (*I* + *T*) and the target signal *v*. The second one (*Jb*) is a background term that consists of the three criteria for 'optimal' mapping given in Equations (2)–(4). These three criteria are used as weak constraints.

$$\begin{array}{lll} J(T) &=& J\_o(T) + J\_b(T) \\ J\_o(T) &=& \|\upsilon - \iota \circ (I + T)\| \\ J\_b(T) &=& \mathcal{C}\_1 \|T\| + \mathcal{C}\_2 \|\nabla T\| + \mathcal{C}\_3 \|\nabla \cdot T\| \end{array} \tag{9}$$

where *<sup>C</sup>*1, *<sup>C</sup>*<sup>2</sup> and *<sup>C</sup>*<sup>3</sup> are three coefficients determined empirically and . is the *<sup>L</sup>*2−norm.

The minimization problem is solved iteratively, for *T* defined on increasingly fine grids. The iterative approach has two advantages. It helps reduce the computational cost and avoids the local minima problem (see below).

In our application, the domain *D* is rectangular. It can be represented by different uniform grids. The regular *nx* × *ny* = *n* grid on which *u*, *v* and *u*morph(*λ*) are given is called the *pixel grid Dn*. The mapping function *T* is defined on a set of coarser grids *Di* (*i* = 1, ..., *I*), called *morphing grids*. It is then represented by two gridded arrays (one for *Tx* and one for *Ty*). The grids *Di* are uniform (2*<sup>i</sup>* + <sup>1</sup>) × (2*<sup>i</sup>* + <sup>1</sup>) = *mi* grids (for *<sup>i</sup>* = 1, ..., *<sup>I</sup>*) covering the domain *<sup>D</sup>*. For *<sup>i</sup>* = 1, ..., *<sup>I</sup>*, the mapping *<sup>T</sup>* discretized on *Di* is noted *Ti*.

The signals *u* and *v*, and so the observation term *Jo* of the cost function, are discretized on the pixel grid *Dn*. The background term *Jb* is discretized on the morphing grid *Di*. We use the second order central scheme except at the boundaries where the first order backward or forward schemes are used. We use bilinear interpolation to estimate the value of *u* and *v* on the distorted grid (e.g., *u* ◦ (*I* + *T*)) and to interpolate *T* on the different morphing grids *Di*.

The finest morphing grid *DI* does not need to be the same as the pixel grid *Dn*. On the contrary, it is computationally advantageous when the morphing grid *DI* has a much coarser resolution. When the number of nodes *mi* of the morphing grids is much smaller than the number of nodes *n* of the pixel grid, solving the minimization problem on the set of morphing grids *Di* is less computationally expensive than to solve it for *T* defined on the high resolution pixel grid *Dn*.
