Phuc Ngo

Yukiko Kenmochi

Akihiro Sugimoto

Hugues Talbot

Nicolas Passat

Discrete rigid registration: A local graph-search approach

Keywords: Image registration, discrete rigid transformation, discrete optimization, discrete rigid transformation graph, graph search

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

Introduction

Image processing and analysis applications are involving an increasing amount of images, which result from several total / partial acquisitions of a same scene or different scenes of same semantics, from different viewpoints and / or times. This multiplicity of data related to a same structure of interest has motivated the development, during the last twenty years, of a novel branch of image processing, devoted to image registration [START_REF] Zitová | Image registration methods: A survey[END_REF]. The most frequent application cases of image registration can be found (non-exhaustively) in medical imaging for inter-or intraindividual mapping from multimodal / multitime acquisition [START_REF] Oliveira | Medical image registration: A review[END_REF], in remote sensing for image orthorectification [START_REF] Schowengerdt | Remote Sensing: Models and Methods for Image Processing[END_REF], in computer vision for stereovision or augmented reality.

Basically, registration consists of mapping "at best" a source image I s onto a target image I t . This implies that I s and I t have equal -or at least similar -semantic contents.

In general, registration techniques involve searching over the parameter space T of certain geometric transformations to find the optimal transformation that maps I s onto I t .

In this context, registration is indeed interpreted as an optimization problem. In addition, it is generally modeled as a continuous problem; in other words, T is considered as a subset of R n . In particular, this continuous paradigm induces an infinity of possible transformations, and the chosen optimization schemes have to handle this issue.

Our working hypothesis is that, since the considered images are generally discrete, and more precisely digital (i.e., defined in Z k mainly with k = 2 or 3), a relevant optimization scheme may consist of navigating within the associated (discrete) parameter space of the transformations from Z k to Z k . In particular, even if such parameter space is huge, it may be computed and explored on-the-flight, leading to combinatorial strategies for determining exact locally optimal solutions, or estimating globally optimal ones.

In this article, we initiate the study of this paradigm, by considering the basic case of rigid transformations between 2D images, i.e., images defined on Z 2 . In the case of rigid transformations, that are obtained by composition of translations and rotations, the parameter space has three dimensions: one for each principal direction of Z 2 , and a third for the rotation angle. The overall parameter space, that is a (finite) subdivision of a subset of R 3 , has a polynomial complexity [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF], and thus cannot be computed as a whole, in practice. Nevertheless, it can be locally computed as a graph. Such a computation, involving a local exploration of the graph, can be integrated into a gradient descent scheme, which leads to a local optimum, for solving the image registration issue. In particular, a multi-scale scheme is proposed in order to speed-up the optimization process. This study then constitutes a proof of concept for a general combinatorial way of considering image registration.

The remainder of this article -that is an extended and improved version of the conference papers [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF][START_REF] Kenmochi | Efficient neighbourhood computing for discrete rigid transformation graph search[END_REF] -is organized as follows. Section 2 introduces and formalizes the optimization problem of image registration, and our specific purpose, namely discrete image registration. Section 3 describes the parameter space of rigid transformations and the way it can be expressed as a combinatorial structure (named a DRT graph [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF]) when handling discrete images. The main contributions of this article can be found in Sections 4 and 5, in which we describe a discrete exploration procedure using neighbouring relations defined on DRT graphs to solve rigid registration issues. Experiments and a conclusion are proposed in Sections 6 and 7, respectively.

Image registration

In general, registration consists of estimating the optimal spatial transformation so that the source and target images are aligned. In other words, it aims to solve an optimization problem expressed as minimizing a metric error between two images, i.e.

T * = arg min T ∈T d((I s • T ), I t ) (1)
where the objective function d to minimize is the distance (dissimilarity) between images. The transformation T * is then searched among the space T of authorised transformations T , defined via the degrees of freedom on their parameters.

Over the years, many approaches have been proposed to register images, each of which is usually designed for specific applications and types of data. All of them practically rely on the following three key-points.

1. "What": the distance between images quantifies how much these images are similar, from a spatial / semantic point of view. In the literature, several distances have been considered for image registration. Depending on the density [START_REF] Hou | Robust dense registration of partial nonrigid shapes[END_REF] / sparseness of images, on the one hand, and intensity-based or feature-based [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF] paradigms, on the other hand, many distances have been introduced, e.g., crosscorrelation, mutual information, ratio image uniformity, least square difference, signed distance, etc. Detailed discussions of these distances can be found, e.g., in [START_REF] Zitová | Image registration methods: A survey[END_REF][START_REF] Boykov | An integral solution to surface evolution PDEs via geo-cuts[END_REF].

"

Where": the geometric transformations that are assumed valid to register the source image onto the target image. The most common transformations are rigid, affine, projective, perspective and global. The determination of these transformations is generally guided by the application and the content of the images.

The range is wide, from rigid to non-rigid transformations; the more complex the transformation, the higher the associated parameter space T (leading to systems that involve from 2 or 3 to thousands of freedom degrees). In addition to geometrical issues, the most complex transformations also induce topological difficulties [START_REF] Noblet | 3-D deformable image registration: A topology preservation scheme based on hierarchical deformation models and interval analysis optimization[END_REF].

3. "How": the optimization scheme provides the modus operandi for estimating the optimal distance by exploring the parameter space of geometric transformations.

The choice of an optimization technique relies, in particular, on the nature of the distance function and the parameter space of transformations. Various continuous optimization schemes have been employed, such as gradient descent [START_REF] Klein | Evaluation of optimization methods for nonrigid medical image registration using mutual information and B-splines[END_REF], conjugate gradient descent [START_REF] Tustison | Directly manipulated free-form deformation image registration[END_REF], Newton-type methods [START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation[END_REF], etc. Due to the ill-condition of the problem, most of these approaches use iterative models of optimization, which consist of performing local search in the parameter space of rigid transformations to find the minima. However, such approaches limit their use in two ways: (1) they strongly rely on the computation of the gradient of the optimizing function, thus this function must be differentiable; (2) they are sensitive to the initial solution and tend to reach local minima. Furthermore, since the problem of image registration is considered in a continuous space, the points of images are assumed to be real values -despite the fact that images are defined in a discrete space. In addition, since the transformation space T is continuous, all possible rigid transformations cannot be reached. More precisely, a sampling process is necessary to explore the search space, potentially biasing to the final result. In particular, the more densely the solution space is sampled, the better accuracy we may achieve, at expense of a higher computational cost.

Nonetheless, there is no guarantee that the obtained result is exact with respect to Equation [START_REF] Zitová | Image registration methods: A survey[END_REF]. Discrete schemes (e.g., graph-cut methods [START_REF] Tang | Non-rigid image registration using graph-cuts[END_REF], particle swarm optimizations [START_REF] Talbi | Particle swarm optimization for image registration[END_REF]) rely on alternative paradigms, where the search space is huge, but finite. Still, the question raised when using these discrete optimizations is how to explore such search space.

We mainly consider this specific issue of discrete schemes for registration, directly related to the third item of the above list. In this context, we focus on the case of 2D images and rigid transformations, i.e., we fix some hypotheses related to the second item. The case of rigid transformations constitutes the simplest -yet non-trivial -case of image registration. This will allow us to emphasise the basics of our combinatorial approach, and provide a sound basis that may be further extended to higher dimension images and transformations. Finally, concerning the distance d for the dissimilarity between images, the choice is -in theory -independent from our optimization framework. Nevertheless, the time cost of distance computation and update has a direct impact on the time cost of the whole process. As a consequence, we will rely preferentially on distances whose local update can be carried out in constant time. The signed distance transform [START_REF] Boykov | An integral solution to surface evolution PDEs via geo-cuts[END_REF], involved in the experiments of Section 6, is an example of such distance.

In our recent works, a fully discrete framework [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF] was proposed to model the space T of rigid transformations over the discrete 2D images. It was proved that the 3D parameter space of the induced discrete rigid transformations is isomorphic to a com-binatorial structure, namely a graph. Navigating within this -so-called Discrete Rigid Transformation (DRT) -graph then allows us to explicitly explore all the potential transformations, by iteratively considering "neighbouring" transformations, that differ pairwise by exactly one point of Z 2 .

In the next sections, we show how such navigation can be carried out in a standard gradient descent paradigm with a low computational cost, by computing on-the-flight the useful part of the parameter space, without building its whole -polynomial sizestructure, thus leading to efficient discrete registration strategies.

Rigid transformations on Z 2 and DRT graph

In order to make this article self-contained, we recall in this section some basic notions of [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF], dealing with rigid transformations defined on R 2 and Z 2 (Section 3.1), the induced parameter space (Section 3.2) and the way this space can be modeled by a combinatorial model (Section 3.3).

Digital rigid transformations

A 2D rigid transformation is defined as a rotation composed with a translation.

In the continuous framework, such a transformation can be formally expressed as a bijective function T : R 2 → R 2 such that for any x = (x, y) ∈ R 2 , the transformed point T (x) has the form

T (x) =           cos θ -sin θ sin θ cos θ                     x y           +           a 1 a 2           (2)
where the parameters a 1 , a 2 ∈ R represent the translation, and θ ∈ [0, 2π[ is the rotation angle. In particular, such a transformation is unambiguously modeled by the triplet of parameters (a 1 , a 2 , θ), and will be often denoted by T a 1 a 2 θ . When applied to an image

I : R 2 → V, it provides a new transformed image I • T : R 2 → V, where V is a finite set of value space. (As examples, if | V | = 2, we say that I is a binary image; if | V | is
equipped with a total order, we say that I is a grey-level image.)

It is not possible to apply directly T to a digital image I : Z 2 → V, since there is no guarantee that T (x) ∈ Z 2 for x ∈ Z 2 . In the discrete framework, the handling of digital rigid transformations requires to define a function T a 1 a 2 θ : Z 2 → Z 2 , which is a discrete analogue of T a 1 a 2 θ . Following a digitization mapping D that associates to each real point x ∈ R its round value [x] ∈ Z, a digital rigid transformation T associated to

T can be conveniently performed by setting T = D • T , as illustrated in the following diagram.

Z 2 T =D•T ------→ Z 2      Id      D R 2 T ------→ R 2 (3)
The function T : Z 2 → Z 2 is then explicitly defined for p = (p, q) ∈ Z 2 by

T (p) = D • T (p) =           [p cos θ -q sin θ + a 1 ] [p sin θ + q cos θ + a 2 ]           (4)
In general, this function is not bijective. However, by setting T -1 :

Z 2 → Z 2 as T -1 = D • T -1 , i.e.
, by considering the standard backward mapping, it is possible to define the digital transformed image I • T -1 : Z 2 → V with respect to T . In the sequel of this article, we focus on such digital rigid transformations. From this point on -for the sake of readability and without loss of correctness -we will note T instead of T -1 , due to the bijectivity of T and T -1 .

From a theoretical point of view, the above notions (images, rigid transformations) are defined on both Z 2 and R 2 . Practically, our purpose is however to study rigid transformations on images of finite size. Under this hypothesis, only some digital rigid transformations are relevant, namely those that actually have an effect on such finite images. We focus on this finite case, and we assume that digital images are defined on subsets of Z 2 of size N × N. Without loss of generality, a digital image I is then written

as I : S → V for S = [[0, N -1]] 2 ⊂ Z 2 .

Digital rigid transformation space

By contrast with rigid transformations defined in R 2 (Equation (2)), their digital analogues are not continuously defined with respect to the parameters a 1 , a 2 and θ. This is pointed out by the fact that a same discrete point T (p) can be obtained with slight variations of (a 1 , a 2 , θ), due to the rounding function D involved in T (Equation ( 4)). In particular, it is shown in [START_REF] Ngo | Topology-preserving conditions for 2D digital images under rigid transformations[END_REF] that digital rigid transformations are generally not bi-150 jective and neither preserve the geometric and topological properties of digital objects.

The above-mentioned discontinuities result in a subdivision of the parameter space The transformations T abθ that lead to the discontinuities of T are expressed, for any p = (p, q) ∈ S that is mapped onto a half-grid point which can be either a vertical A set of tipping surfaces induces a subdivision of the (a 1 , a 2 , θ) space into equivalence classes, each of which consists of rigid transformations leading to a same digital result in Z 2 . Such classes are called discrete rigid transformations (DRTs).

R 3 of (
p Φ = (k + 1 2 , λ) ∈ H, or a horizontal p Ψ = (λ, l + 1 2 ) ∈ H (with k, l ∈ Z and λ ∈ R) by Φ pqk : R 2 -→ R (a 2 , θ) -→ a 1 = φ pqk (θ) = k + 1 2 + q sin θ -p cos θ (5)
Ψ pql : R 2 -→ R (a 1 , θ) -→ a 2 = ψ pql (θ) = l + 1 2 -p sin θ -q cos θ (6)

Discrete rigid transformation graph

By mapping each 3D open cell onto a 0D vertex, and each separating 2D surface onto an 1D edge -following a Voronoi / Delaunay duality -we can model the subdivision of the 3D parameter space of (a 1 , a 2 , θ), as a graph, called DRT graph.

Definition 1 (DRT graph [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF]). Let G = (V, E) be the graph defined such that:

(i) any vertex v ∈ V models a 3D open cell associated to a DRT; (ii) any edge e = (v, w) ∈ E connects two vertices v, w ∈ V sharing a 2D surface. The graph G = (V, E) is called a DRT graph.
In other words, the parameter space R 3 of rigid transformations is subdivided, and

the set V of G = (V, E)

Discrete rigid transformation graph for image registration

DRT formulation of rigid image registration

On the one hand, the registration problem stated in Section 2 and formalized in Equation ( 1) consists of finding the transformation that allows us to map at best a source image I s onto a target image I t with respect to a given distance. On the other hand, the DRT graph G = (V, E) associated to an image of finite size, models all the existing DRTs for this image, i.e., all its transformed images; in particular this set of solutions is finite.

Therefore, one can use G to find the best rigid registration between I s and I t . More precisely, let G = (V, E) be the DRT graph associated to the source image I s . There

exists v ∈ V such that d(I s • T v , I t ) is minimal, with T v the transformation associated
to the vertex v and I s • T v the transformed image of I s associated to T v 5 . By that way, Equation ( 1) can be reformulated by considering the rigid transformation space

T = {T v | v ∈ V}
of all the DRTs as follows:

T * v = arg min v∈V d(I s • T v , I t ) (7)
Roughly speaking, a brute-force search within the whole graph G would lead to a global optimal solution for Equation [START_REF] Hou | Robust dense registration of partial nonrigid shapes[END_REF]. In particular, from Section 3.3, the DRT graph transforms progressively (i.e., pixel by pixel) the image. Therefore, the exhaustive scanning of the DRT graph can be done by updating incrementally the distance between images.

However, it is proved in [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF] that the DRT graph G, associated to an image of size N × N, has a space complexity of O(N 9 ), which induces the same time complexity both for construction and exhaustive exploration of G. This polynomial complexity practically forbids the use of G as a whole in real applications for large images.

To overcome this issue, a local approach is proposed to reach a local optimum in G.

Such an approach is based on a local search allowing to efficiently explore the graph without generating it entirely, and thus leads to a much lower algorithmic complexity.

Thanks to the property of local modifications along the edges of the DRT graph, this approach can be integrated into an optimization process that guarantees to reach such optimum in the sense of minimizing the objective function with respect to a given distance between images.

Local search in a DRT graph

In this section, we propose a combinatorial optimization -gradient descent -applied on the investigated structure of DRT graph, in order to find locally optimal solutions for the rigid registration problem stated above. To this end, we first need to introduce a notion of neighbour of DRTs, and properties related to the locally optimal solutions induced by such neighbourhoods.

Neighbourhood in a DRT graph

In order to perform a gradient descent, for searching a local optimum, two questions are raised: (i) the choice of the initial vertex / transformation; and (ii) the choice of a "good" search area around a vertex, namely neighbours of this vertex. The latter is the most critical issue while the first is often guided by the application.

From the above discussions on the definition, structure and semantics of a DRT graph, we define the notion of neighbours by using the standard definition in graph theory. Two DRTs are adjacent -neighbours -in a DRT graph if their respective vertices are adjacent (i.e., share an edge) in this graph. Equivalently, the two DRTs share a surface in the dual parameter space of (a 1 , a 2 , θ). Based on the adjacency relation in the DRT graph, we can define the neighbours of depth k ≥ 1 of a vertex, namely k-neighbours, as follows.

Definition 2 (DRT neighbourhood). Given a DRT graph G = (V, E), the k-neighbour of a vertex v ∈ V is N k (v) = N k-1 (v) ∪ u∈N k-1 (v) {w ∈ V | ∃e = (u, w) ∈ E} where N 0 (v) = {v}.
This notion is illustrated in Figure 3.

Gradient descent within a DRT graph

We now describe a generic local search algorithm, namely a gradient descent, Output: A set V ⊆ V of local optima with respect to Equation [START_REF] Hou | Robust dense registration of partial nonrigid shapes[END_REF].

1 V, V * ← {v 0 } 2 while ( V = V * ) do 3 V ← V * 4 V * ← arg min d(I s • T v , I t ) | v ∈ v ′ ∈V * N k (v ′ )
we can then find the vertices presenting the minimal value d(I s • T v , I t ). The obtained DRT(s) then become(s) new seed(s) for the next step.

By iterating this process, it is then possible to carry out a gradient descent from a seed vertex v 0 ∈ V, that finally leads to a local optimum with respect to Definition 2.

250

The number of vertices of the DRT graph G being finite, this process, summarized in Algorithm 1, has a guaranteed termination.

Due to digitization effects, the distance d may be non-smooth, and can locally increase on a globally decreasing trajectory. In order to avoid the side effects on such disturbance on the gradient descent, that requires a globally decreasing distance along a reasonably long trajectory, it is preferable to compute ∆d on a k-neighbourhood, with k > 1, to improve the robustness of the process.

Computing the 1-neighbourhood N(v) of a vertex v in low time cost is indeed tractable in a DRT graph, but preserving a good time complexity for larger neighbourhoods N k (v) constitutes a real challenge. In particular, it is the main issue to tackle in our proposed framework.

In [START_REF] Ngo | Combinatorial properties of 2D discrete rigid transformations under pixel-invariance constraints[END_REF], we initially proposed an algorithm which could be adapted to find 1neighbourhoods N(v). This algorithm presents a complexity of O(N 4 log N), where N × N is the size of the given image. Based on the relations that link tipping surfaces and tipping curves (Section 3.2), a better algorithm was proposed in [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF] to compute

N(v) with a complexity of O(mN 2 )
, where m is the maximum degree of the DRT graph (i.e., the number of 1-neighbours of a DRT). A recursive version of this algorithm could be derived for computing k-neighbourhoods, but it has an exponential complexity O(m k N 2 ) with respect to k. In the next section, we present a more efficient algorithm, initially introduced in [START_REF] Kenmochi | Efficient neighbourhood computing for discrete rigid transformation graph search[END_REF], that allows us to compute N k (v), with a low complexity

O(kN 2 ).
It should be mentioned that the proposed framework of registration based on DRT graph search is valid for any digital 2D image presenting isotropic pixels, since DRT graphs, by definition, do not depend on the values that are assigned to the pixels of images. In other words, the structure is invariant for any image defined on the same support; i.e., any image of same size N×N (see Section 3.3). Furthermore, the proposed framework is also independent from the choice of a distance metric, which is calculated at each iteration of the gradient process. For the sake of readability -but without loss of generality -we will consider binary and grey-level images with the distance based on the signed distance function [START_REF] Boykov | An integral solution to surface evolution PDEs via geo-cuts[END_REF][START_REF] Molchanov | Distance transforms for real-valued functions[END_REF] in the experiences of Sections 6.3 and 6.5.

Neighbourhood construction

The complexity of the above local search algorithm (Algorithm 1) directly relies on the complexity of the construction of the k-neighbours N k (v), i.e., the computation of a local part of the DRT graph modeling the neighours of a given vertex v. In this section, we describe an efficient algorithm for constructing such neighbours, based on the sweeping plane technique introduced in [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF]. To this end, we first describe a modified sweeping algorithm for DRT sub-graph construction in an interval of θ within [0, 2π[.

Sweeping algorithm for discrete rigid transformation sub-graph construction

Similarly to the algorithm proposed in [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF] for building incrementally the whole DRT graph, we describe hereafter how to build partially the graph within a given range of θ values, smaller than [0, 2π[. This problem is formalized as follows: given a set S of In the sequel, we describe how to use this partial sweeping technique for computing the direct neighbours N k (v) of v, of depth k = 1 and then k ≥ 1. To guarantee the efficienty of this method, we will need as inputs: (1) the analytic representation of a DRT, and (2) a "representative" continuous rigid transformation within a given DRT.

Analytic representation of discrete rigid transformations

Let I v be the transformed image associated to the vertex v ∈ V, obtained by applying the rigid transformation(s) of v on the initial image I. Each pixel (p i , q i ) of I v has a corresponding pixel (p ′ i , q ′ i ) in I. Such correspondence is modeled by the following inequalities deriving from Equation ( 4)

p ′ i - 1 2 < p i cos θ -q i sin θ + a 1 < p ′ i + 1 2 ( 8
)
q ′ i - 1 2 < p i sin θ + q i cos θ + a 2 < q ′ i + 1 2 (9)
Since I has a size N × N, each of its N 2 pixels generates 4 such inequalities. Then, this set of 4N 2 inequalities provides an analytic expression of the DRT v. In the parameter space (a 1 , a 2 , θ), the obtained 4N 2 inequalities then define a 3D cell, called feasible region and denoted as R v . In other words, such region gathers all the parameter triplets associated to the DRT corresponding to the vertex v.

When interpreting these inequalities in terms of tipping surfaces/curves (Equations (5-6)), it appears that for each pixel of I v , Equations (8-9) define a region of the parameter space that is bounded by two offset vertical (resp. horizontal) tipping surfaces/curves Φ p i q i p ′ i and Φ p i q i p ′ i -1 (resp. Ψ p i q i q ′ i and Ψ p i q i q ′ i -1 ). For any i ∈ [[1, N 2 ]], Φ p i q i p ′ i (resp. Ψ p i q i q ′ i ) is called an upper tipping surface/curve, while Φ p i q i p ′ i -1 (resp.

Ψ p i q i q ′ i -1
) is called a lower tipping surface/curve. The sets composed by these surfaces/curves, for all i ∈ [[1, N 2 ]], are denoted S + 1 (I v ) and S - 1 (I v ) (resp. S + 2 (I v ) and S - 2 (I v )), and exemplified in Figure 5.

Representative of a discrete rigid transformation

From the definition of DRTs (Section 3.2), the feasible region R v of a DRT v contains the set of all the rigid transformations leading to a same digital transformation.

We call a representative of R v one of these rigid transformations, further noted T v .

From the sweeping algorithm described above, whenever the sweeping plane γ reaches an intersection in the parameter space, a 3D cell is ended while a new is generated. Equivalently in the dual space, new vertices and edges are generated in the DRT sub-graph. By tracking these evolutions during the sweep, one can determine the θ min

S 2 = {S + 2 (I v ) ∪ S - 2 (I v )} (Section 5.
2), the θ min (resp. θ max ) of the feasible region R v of v is determined at the θ value for which the sweeping plane starts (resp. finishes) intersecting the tipping curves of each plane (a * , θ) in the sequences of upper and lower curves (Figure 6).

In [START_REF] Ngo | On 2D constrained discrete rigid transformations[END_REF], it was observed that the feasible regions of DRTs present directional convexity properties with respect to the a ⋆ -axes in the parameter space (a 1 , a 2 , θ). Thanks to these properties, we can calculate a representative once knowing its θ min and θ max6 as follows.

T v = (a 1v , a 2v , θ v ) associated to R v
                 θ v = θ min (v) + θ max (v) /2 a 1v = arg min { φ + (θ v ) | φ + ∈ S + 1 (I v ) } + arg max { φ -(θ v ) | φ -∈ S - 1 (I v ) } /2 a 2v = arg min { ψ + (θ v ) | ψ + ∈ S + 2 (I v ) } + arg max { ψ -(θ v ) | ψ -∈ S - 2 (I v ) } /2 (10)
As stated above, T v can be any transformation within R v . Thus, θ v can be any value in

the interval ] θ min (v), θ max (v) [ of R v .
In particular, we can set θ v at the middle of the interval as in Equation [START_REF] Noblet | 3-D deformable image registration: A topology preservation scheme based on hierarchical deformation models and interval analysis optimization[END_REF].

In [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF], it is shown that the value of θ can be expressed and compared as a periodic continued fraction, modelled by a sequence of integers. Using this representation, one can find a θ v between θ min (v) and θ max (v) in Equation 10in an exact way.

It may be noticed that in the sweeping algorithm for DRT sub-graph construction (Section 5.1), we only need the interval of θ values as parameter, namely θ start and θ end (with θ start = θ v and θ end = either θ min or θ max ). In such cases, the explicit calculation of a 1v and a 2v is not required.

1-Neighbours of a discrete rigid transformation

From Definition 2 and Figure 3, computing the 1-neighbours of a DRT v is equivalent to finding the vertices adjacent to v. Let us assume that we know a representative T v associated to v (calculated as shown in Section 5.3, or given as input). We can now formulate the problem as follows: given a rigid transformation T v associated to the

feasible region R v of a DRT v, represented by 4N 2 inequalities of S 1 = {S + 1 (I v ) ∪ S - 1 (I v )} and S 2 = {S + 2 (I v ) ∪ S - 2 (I v )} (Section 5.
2), how can we determine the neighbourhood N(v) of v.

The computation of N(v) is equivalent to constructing the DRT sub-graph of only neighbouring vertices and edges of v. From Section 5.1, we already have a sweeping algorithm for building the DRT sub-graph for θ ∈ [θ start , θ end [. We now use it for finding N(v), by sweeping the plane γ twice. Indeed, the value7 θ v is considered as the initial position, and the plane is swept to the left-hand and right-hand sides of θ v until the θ min and θ max values of R v are reached, respectively (Figure 7). At each step of the algorithm -when the sweeping plane reaches an intersection of tipping surfaces -several vertices and edges are generated, but only those that are adjacent to v are considered for the DRT sub-graph of N(v).

The question raised by using the algorithm of Section 5.1 is how to limit the construction of the DRT sub-graph to these neighbours only. In fact, thanks to the directional convexity property along the a ⋆ -axis of a feasible region of a DRT, the sequence of tipping surfaces intersecting the sweeping plane is separated into two parts: the upper and the lower tipping surfaces, respectively (Section 5.2). Thus, by observing this sequence of tipping surfaces updated during the sweep, one can track them until the algorithm stops. The procedure for the left-hand side is given in Algorithm 2 (the right-hand side is similar).

k-Neighbours construction algorithm

In this section, we extend the previous algorithm for finding the k-neighbours of a given DRT v, i.e., N k (v) for k ≥ 1 (Figure 8). For this, we need a function that associates the depth of a vertex u with respect to v, denoted by t v (u) (in particular,

t v (v) = 0).
For each new vertex u generated during the sweep, the depth of u is defined from its two adjacent vertices w 1 and w 2 in the a ⋆ -direction of the tipping surface intersection, such that t v (u) = 1 + min{t v (w 1 ), t v (w 2 )}. (Note that an iterative backtracking process Algorithm 2: 1-Neighbours computation (left-hand side along the θ-axis)

Input: A DRT v / its associated image I v .

Output:

The DRT sub-graph G v = (V v , E v ) containing the 1-neighbours of v. 1 Initialize V v = ∅, E v = ∅ 2 Initialize θ prev = θ v 3 for ⋆ = 1, 2 do 4
Determine the tipping surfaces S + ⋆ (I v ), S - ⋆ (I v ) associated to v (Section 5.2).

θ next = min{θ + 1 , θ - 1 , θ + 2 , θ - 2 } ∆G v = S weep(S 1 ∪ S 2 , θ prev , θ next ) if ∃u ∈ ∆V v such that e = (u, v) ∈ ∆E v then G v = G v ∪ ∆G v θ prev = θ next if the next intersecting surface with f + ⋆ (or f - ⋆ ) is in S ⋆ then
Exchange their order in S ⋆ .

else

Replace f + ⋆ (or f - ⋆ ) in S ⋆ with the new intersecting surface.
until no more separation of tipping surfaces in S ⋆ ; is needed to update the depths of w ⋆ and its successive neighbours, whenever t v (w ⋆ ) >

t v (u) + 1.)
This depth function is used to terminate the procedure of construction of the DRT 405 sub-graph N k (v). Indeed, if t v (u) > k for all vertices in the current sweeping plane, then the process ends. Furthermore, in order to obtain all the vertices u such that t v (u) ≤ k, we need to initialize the set of tipping surfaces S + ⋆ (I v ) (resp. S - ⋆ (I v )) with the k + 1 lowest (resp. uppermost) tipping surfaces that intersect the initial sweeping plane at θ v .

Note that in the DRT sub-graph of N k (v), we keep only vertices whose neighbourhood depth is lower than k. The modified algorithm for N k (v) is given in Algorithm 3. This algorithm is similar to Algorithm 2, except for three points (depicted in blue):

(1) the initialization of the set S ⋆ with 2(k + 1) tipping surfaces instead of 4 (Lines 3-6); (2) the determination of vertices whose neighbourhood depth is ≤ k, to add them into the sub-graph of N(v) (Line 12); and (3) the second loop termination when all considered vertices have a depth > k (Line 19). Roughly speaking, Algorithm 2 is a particular case of Algorithm 3 for k = 1.

Algorithm 3: k-Neighbours computation (in the left-hand side along the θ-axis)

Input: A DRT v / its associated image I v ; a positive integer k.

Output:

The DRT sub-graph G v = (V v , E v ) containing the k-neighbours of v. 1 Initialize V = ∅, E = ∅ 2 Initialize θ prev = θ v 3 for ⋆ = 1, 2 do 4
Determine the tipping surfaces S + ⋆ (I v ), S - ⋆ (I v ) associated to v (Section 5.2).

5

In S + ⋆ (I v ) (resp. S - ⋆ (I v )), find the (k + 1)-th lowest (resp. uppermost) tipping surface f + ⋆ (resp. f - ⋆ ), that intersects the initial plane at θ v .

6

Find and sort the k + 1 tipping surfaces that are upper (resp. lower) or equal to f + ⋆ (resp. f - ⋆ ), and put them in an ordered set S ⋆ .

In Algorithm 3, the first loop (Lines 3-6) initializes the set S ⋆ used to compute the k-neighbours N k (v) of a given DRT v. More precisely, Line 4 finds the sets S + ⋆ (I v ) and S - ⋆ (I v ) containing the upper and lower tipping surfaces of v, and has a complexity of O(N 2 ). Lines 5-6 compute the 2k nearest tipping surfaces after the lowest and uppermost tipping surfaces f + ⋆ and f - ⋆ and store all in S ⋆ . (Note that S ⋆ contains 2(k+1) tipping surfaces.) In particular, O(N 2 ) for Line 5 in average, if we use Hoare's find algorithm [START_REF] Hoare | Algorithm 65: Find[END_REF]; and O(N 2 ) and O(k log k) for finding and sorting the tipping surfaces in Line 6, respectively.

The second loop (Lines 7-18) calculates N k (v) from S ⋆ . Firstly, Lines 8-10 compute the interval of interest of θ, by calculating the next intersection (namely, θ next ) from the current θ (namely, θ prev ). In this step, the most costly part is Line 9, which finds all intersections of S ⋆ . Since S ⋆ contains 2(k + 1) tipping surfaces, this costs

O(|S ⋆ | 2 ) = O(k 2 )
. Secondly, Line 11 builds a DRT sub-graph by using the Sweep algorithm for the interval [θ prev , θ next [; this requires O(N 2 ), which is the most costly part in the loop. The number of iterations of this second loop is estimated as m(2k + 1), since the k-neighbours of the vertex v contain at most 2k + 1 adjacent vertices in the θ-direction. Therefore, the complexity of Algorithm 3 is O(mkN 2 ). One remarks that this complexity depends on the size of the neibourhood of v, namely the value of m. We now prove that m has a low constant value, that leads to a final complexity of O(kN 2 ).

From [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF], the DRT graph space complexity for an image of size N × N is O(N 9 ), both for vertices and edges, since the number of vertices and that of edges grow at the same rate. However, we can infer that m is actually bounded, independently of N. Indeed, by analogy, let us imagine that we divide a 2D plane with straight lines defined randomly. Three lines will almost never intersect at a same point, and for a number of lines sufficiently large, the cells of the induced subdivision will be mostly triangles. Following this analogy, we may infer that the degree of the vertices of the 2D DRT graphs in the planes (a 1 , θ) and (a 2 , θ) is close to 3, in average. However, this analogy has some limits: the considered tipping curves are not straight lines, while their very regular structure implies that many curves often intersect at a same point.

Nevertheless, we can assimilate a 2D DRT graph (which is the projection of a 3D DRT graph onto the (a ⋆ , θ) plane) to a planar graph whenever N is sufficiently large. Under such assumption, the Euler formula is valid, i.e., we have ve + f = 2, where v, e and f are the number of (0D) vertices, (1D) edges and induced (2D) cells, respectively.

From the construction of the DRT graph, each cell f is bounded by at least 4 edges e 455 [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF]. In other words, we have 4 f ≤ 2e. It then comes that 2e/v ≤ 4 -8/v. As v ≫ 8 in DRT graphs, we have 2e/v < 4, where 2e/v is indeed the average degree of the 2D DRT graph. It follows that the average degree m of the 3D DRT graph (obtained by Cartesian product of two 2D DRT graphs) is lower than 2 × 4 = 8. Figure 9 illustrates the distribution of these degrees for a 2D DRT graph. 460

Experimental results

In this section, we provide experiments devoted to validate the theoretical complexity analysis described in Section 5.6, as well as to demonstrate the behaviour of the local search approach using the algorithms developed in Section 5, in the context of image rigid registration.

Average degree of DRT graphs

First, we observe the evolution of the average degree of the vertices within a DRT graph, i.e., the number m of 1-neighbours of a DRT. From Section 5.6, we proved that the average value of m is lower than 8. This is confirmed by the experimental analysis illustrated in Figure 10. The value m is calculated using Algorithm 2. The experiments are carried out for images of size varying from 5 ×5 to 80 ×80, in each of which several DRTs are chosen randomly to compute the value of m.

Observing the experiments on 2D DRT graphs -the projections of a 3D DRT graph onto the (a ⋆ , θ) planes -the degree of the vertices remains close to the average value, i.e., 4 (this is 8 in case of 3D DRT graphs, which is twice of the 2D one). Indeed, Figure 11 shows that the 2D DRT vertex degrees converge, with respect to the image sizes, to a stable distribution that contains mainly degrees of value 3 and 4 (with a maximal value experimentally identified at 8). This confirms the stable behaviour of m with respect to the image size, in practice.

Run-time complexity

Experiments: Image registration

We now describe some experiments and results of image registration, using an iterated local search developed in a fully discrete framework. More precisely, we combine a gradient descent method (Section 4.2) for the optimization part of the scheme, with the k-neighbours computation (Section 5.5) for the local search part. For run-time efficiency of Algorithm 1, instead of maintaining a set of all local optima at each gradient descent iteration, we will select randomly one in the set to iterate. It should be noticed that the proposed approach can be carried out on various types of images (binary, greylevel, label, color, etc.) since DRT graphs are defined independently from the value space of the images (Section 3.3). It is however necessary to use an appropriate distance function d in Equation [START_REF] Hou | Robust dense registration of partial nonrigid shapes[END_REF]. In order to illustrate and analyze the issues related to the digitization on the discrete space of the transformed images -and for the sake of readability -we focus our experiments on the case of binary and grey-level images with with the distance based on the signed distance function [START_REF] Boykov | An integral solution to surface evolution PDEs via geo-cuts[END_REF][START_REF] Molchanov | Distance transforms for real-valued functions[END_REF]. The latter choice is justified as follows.

We recall that the signed distance describes the distance of a given point to the boundary of the object. More precisely, the distance is positive (resp., negative) if the point is inside (resp., outside) the boundary, and its absolute value decreases towards this boundary, where the signed distance function is zero (see Figure 13). Formally, let Ω ⊂ Z 2 be a bounded domain, for any point x ∈ Z 2 , the signed distance function

sd f : Z 2 → R is defined as: sd f Ω (x) = d Ω (x) -d Ω (x)
where d Ω (x) denotes the usual Euclidean distance transform to the set Ω and Ω is the complement of Ω. Then, the distance based on signed distance function between two bounded sets Ω and Ω ′ , which are in the target image I t and the transformed source image I s • T respectively, is defined as:

d(Ω, Ω ′ ) = Ω ′ sd f Ω (x) - Ω sd f Ω (x)
In the context of registration, the second term is constant, while the first term is 505 defined as a point-wise distance between the transformed image of the source and the distance map of the target -generated once for the whole process 8 . From Section 3.3, two adjacent/consecutive transformed images in the DRT graph differ by only one pixel. By that way, at each iteration of the gradient descent method, we need to update the distance only at this changed pixel with respect to the distance map of the 510 target image. In other words, the distance can be computed in a constant time, except for the very first iteration when the distance is calculated for the whole image. Furthermore, it is shown in [START_REF] Boykov | An integral solution to surface evolution PDEs via geo-cuts[END_REF][START_REF] Molchanov | Distance transforms for real-valued functions[END_REF] that this distance gives fewer flat zones for the gradient term of d.

Concerning the initial solution, namely a seed, for the local search algorithm, we use the first and second order central moments for binary images [START_REF] Flusser | Moments and Moment Invariants in Pattern Recognition[END_REF], and a SIFT feature based method [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF][START_REF] Amintoosi | A fast image registration approach based on SIFT key-points applied to super-resolution[END_REF] for grey-level images.

The experiments are carried out with different neighbourhood sizes k = 1, 3, 5, 10, 15 on binary and grey-level images of size 53 × 53 from given initial transformations (Figures 14 and16). The results obtained with the proposed algorithm are shown in Figures 15 and17, respectively. We can observe from these experiments that the locally optimal distance is improved when we use a larger neighborhood, which is coherent in a gradient descent paradigm.

Figure 18: Multi-scale approach for discrete rigid transformations. The image is considered in blocks of pixels of size [-3 2 , 3 2 ] 2 , and a resolution of 1 3 . In the parameter space, at resolution 1 3 we obtain the tipping surfaces / curves depicted by continuous lines. When increasing the resolution, "offset" surfaces / curves (dashed lines) appear.

Multi-scale handling of the discrete rigid transformation graph

The gradient descent approach proposed for handling the rigid registration issue requires to compute, at each iteration of the descent, the k-neighbours of the current DRT. The cost of the descent is then O(skN 2 ), where s is the number of iterations until reaching a local minimum. On the one hand, when s is small, a local minimum is reached rapidly, with a high probability to converge on a non-satisfactory result, unless the initial DRT is chosen "close" to the optimal solution. On the other hand, when s is large, the overall computational cost is high.

In order to allow a robust, non-costly gradient descent registration, without requiring an accurate initialisation of the process, a solution consists of considering a multiscale variant of the proposed scheme. This is tractable since the DRT graph presents a structure that is compliant with multiscale decomposition.

More precisely, for a scale parameter 1 3 , i.e., by considering "superpixels" on intervals [- 3 2 , 3 2 ] 2 instead of [-1 2 , 1 2 ] 2 , the tipping surfaces / curves of the initial DRTs are "subsampled" by a factor 3, as illustrated in Figure 18.

Based on this fact, a first gradient descent can be considered on a coarse image of size N/3 n × N/3 n , sampled from the real image, with a computational cost : Distance evolution at each scale level during the multi-scale local search for the inputs in Figure 14 using the same seed as in Figure 14.

O(sk(N/3 n ) 2 ). The local minimum obtained at this scale can then be used as initialization for a new gradient descent in a finer image of size N/3 n-1 × N/3 n-1 with a cost O(sk(N/3 n-1 ) 2 ), and so on, until processing the real image of size N × N, with a cost O(skN 2 ). By that way, for an average number of s iterations at each scale, the overall "distance" on the DRT graph is s. n k=0 3 k = s.O (3 n+1 ), for a computational cost 545 of n.O(skN 2 ). The speed-up of this multiscale approach is then 3 n+1 /n, i.e., exponential compared to the standard gradient descent. In addition, it allows to initialize the process potentially far from the optimal solution.

Figure 19 shows the experimental results using the multi-scale approach for the source and target images in Figure 14

Application case: Registration of granular material images

In [24], a study on the relative displacement of grains, under directional shear, was investigated. To this end, several experiments were performed, in which an assembly of granular objects is used to simulate the grains. The object volume was slowly moved in simple shear at constant stress, as illustrated in Figure 21. Though under a homogeneous deformation, the granular objects do not have a homogeneous displacement.

The objective of the experiments was to observe the behaviour, in terms of translations and rotations, of each granular object under such shear stress. To this end, 2D pictures were taken at different moments of the experiment. These pictures were used to track and calculate the rigid motion from one image to another.

We use the proposed method to perform rigid registration on such granular objects over the images. It should be noticed that segmentation and tracking of the objects were necessary before applying our registration approach. Figures 22 and23 show the results obtained for the input of Figure 21.

Still in [24], a digital image correlation approach was used to determine the motions of each granular object over the experiments. Using this result as initial transformations of our approach, Table 1 shows that the final results can be slightly improved by the proposed approach.

Conclusion

We have proposed a first attempt for a fully discrete framework of image registration under rigid transformations based on DRT graphs. Using a DRT graph, one can potentially generate all the transformed images of a given image under registration.

575

Such a generation of transformed images presents a complexity in linear time with re- This first algorithmic and experimental discussion opens the way to further developments of combinatorial approaches of image registration. On the one hand, the ways to extend this framework of rigid registration from 2D to 3D images will be investigated. This will require, in a first time, to accurately understand the subtleties of the analogue of the DRT graph for 3D transformations, and to solve some analytical open

  (a) Tipping surfaces. (b) Tipping curves.

Figure 1 :

 1 Figure 1: (a) Tipping surfaces in the 3D parameter space (a 1 , a 2 , θ) and (b) their cross-sections, namely tipping curves, in the 2D planes (a 1 , θ) and (a 2 , θ).

  (a) Parameter space. (b) DRT graph.

Figure 2 :

 2 Figure 2: (a) Subdivision of the (a 1 , a 2 , θ) parameter space into 3D cells by tipping surfaces and (b) the associated DRT graph.

  is the direct transposition of this subdivision. In particular, G associated to a digital image I of finite size N × N is a finite data-structure that describes all the possible digital rigid transformations of I. An example of DRT graph is illustrated in Figure2.It should be mentioned that the DRT graph G does not contain any geometric parameters (a 1 , a 2 , θ) of rigid transformations but only the topological information, that models the relationship between any neighbouring transformed images. More precisely, each vertex v ∈ V of G is associated with a unique transformed image generated by any digital rigid transformation in v. It is showed in[START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF] that two connected transformed images in G differ by at most one over the N 2 pixels of I. Thus, each edge e = (v, w) ∈ E of G is -implicitly -associated to a function indicating the only modification that differs between the transformed images corresponding to the DRTs v and w. This allows us to use the DRT graph to produce all the transformed images via successive elementary (i.e., single-pixel) modifications. This property opens a way of involving the DRT graph in digital image processing and analysis tasks. In the next section, we show how to use this combinatorial structure for handling image registration.

  within the DRT graph by exploring k-neighbourhoods of DRTs. Given two images I s and I t , the DRT graph G = (V, E) associated to I s , and an initial DRT v ∈ V (namely a seed), we determine a local optimum (or a set of local optima), i.e., a DRT minimizing the distance d between I s and I t . Since each DRT corresponds to a unique transformed image, and as there is at most one different pixel between two neighbouring DRTs (Section 3.3), it is possible to locally and incrementally compute the gradient ∆d of d between the transformed image obtained from I s and the target image I t . By examining ∆d for all neighbouring vertices of N(v) for a given vertex v, (a) Parameter space. (b) Neighbours in the DRT graph.

Figure 3 :Algorithm 1 :

 31 Figure 3: Example of k-neighbours of a DRT. The given DRT is depicted in red, its 1-neighbours in green and its 2-neighbours in green and yellow.

s 1 Figure 4 :

 14 Figure 4: DRT graph construction by the sweeping plane algorithm, with 2 vertical (blue, cyan) and 2 horizontal (red, magenta) tipping surfaces. (a) 3 × 3 rectangular cells generated by the tipping surfaces in each sweeping plane. (b) The associated part of the DRT graph in each plane (in green: new vertices and edges in the second and third planes).

Figure 5 :

 5 Figure 5: Analytic representation of a DRT by its set of tipping surfaces. The figure illustrates the projection of these surfaces onto the planes (a ⋆ , θ), where the upper (resp. lower) tipping surfaces/curves are represented by continuous (resp. dashed) lines, while the feasible region of the DRT is depicted in green.

Figure 6 :

 6 Figure 6: A representative T v = (a 1v , a 2v , θ v ) of the feasible region(s) R v (in green) associated to a DRT v, is any rigid transformation within R v . Note that any sweeping plane (in grey) orthogonal to the θ-axis, and passing through R v intersects the tipping curves of each plane (a 1 , θ) and (a 2 , θ) in a sequence in which the upper and lower curves are separated.

Figure 7 :

 7 Figure 7: Illustration of 1-neighbours construction of a given DRT v on the plane (a 1 , θ). The region of v is in green while its 1-neighbours are in red.

Figure 8 :

 8 Figure 8: Illustration of k-neighbours (k = 2) construction of a given DRT v on the plane (a 1 , θ). The region of v is in green, while its 1-and 2-neighbours are in red and purple, respectively.

Figure 9 :

 9 Figure 9: Degree distribution in a 2D DRT graph, viewed in the dual subdivision of the parameter space.Each colour represents a given degree, that corresponds here to the number of curves bounding each cell (3: red, 4: green, 5: blue; 6: yellow).

Figure 10 :

 10 Figure 10: Degree of vertices m within DRT graphs. The value m is calculated over 960 experiments for different image sizes. Red (resp. blue and green) points denote the average value of m (resp. the standard deviation ±σ) over the experiments.

Figure 11 :

 11 Figure 11: Normalized vertex degree distribution in 2D DRT graphs.

480

  We now observe the run-time computation of the k-neighbours algorithm. As described in Sections 5.5 and 5.6, Algorithm 3 has a theoretical complexity of O(kN 2 ), which is in linear time with respect to k for a given image size N × N, and in linear time with respect to this size N × N, for a given k. We experimentally assess the actual cost of this algorithm. To this end, we firstly consider an image of size 53 × 53 and 485 vary the value of k from 1 to 15. Secondly, we fix k-neighbour size at 3 and vary N from 10 to 50. The results are shown in Figure 12; for each varying value of k and N, the experiments are repeatedly performed, and the time values are then averaged, to confirm such linear result.

Figure 12 :

 12 Figure 12: Run-time complexity of k-neighbours computation is experimentally in linear time with respect to k (left) and image size N 2 (right). Red (resp. blue and green) points denote the average value (resp. the standard deviation ±σ) of time over the experiments.

Figure 13 :

 13 Figure 13: Examples of signed distance map of images. (a) Binary image. (b) Signed distance transform of (a). (c) Grey-level image. (d) Average of the signed distance function of all binary level set images of (c).

  image (c) Initial solution (d) Solution: k = 1 (e) k = 3 (f) k = 5 (g) k = 10 (h) k = 15

Figure 14 :

 14 Figure 14: Input images and results of the iterated local search for image registration. (a) Source image, (b) Target image generated from (a) with (a 1 , a 2 , θ) = (0.1, 0.3, 0.114) and (c) the initial transformed image of (a) with (a 1 , a 2 , θ) = (0.365, -0.045, 0.1423). (d-h) Local optima obtained from (c) by using k-neighbours for k = 1, 3, 5, 10, 15; the representative transformation parameter triples of the resulting DRT are (a 1 , a 2 , θ) = (0.3641, -0.0483, 0.1423), (0.0118, 0.0176, 0.1458), (-0.0109, 0.0125, 0.0147), (0, -0.0529, 0.1446), (-0.0385, -0.0619, 0.1442) respectively. Note that in (c-h), pixels are coloured if they are different from those in (b); yellow (resp. red) pixels are considered as white (resp. black) in (c-h) and black (resp. white) in (b).

Figure 15 :

 15 Figure 15: Distance evolution during iterations of local search for the inputs in Figure 14 (a) and (b), from the initial transformation in Figure 14 (c), with respect to different depths k.

Figure 16 :

 16 Figure 16: Input images and results of the iterated local search for image registration. (a) Source image, (b) Target image generated from (a) with (a 1 , a 2 , θ) = (0.1, 0.3, 0.114) and (c) the initial transformed image of (a) with (a 1 , a 2 , θ) = (0.49, 0.52, 0.153) and (d) difference between the target image (b) and the initial transformed image (c). (e,g,i,k) Local optima obtained from (c) by using k-neighbours for k = 1, 5, 10, 15; the representative transformation parameter triples of the resulting DRT are (a 1 , a 2 , θ) = (0.491, 0.5204, 0.153), (0.4502, 0.5143, 0.1531), (0.4529, 0.5118, 0.1535), (0.4673, 0.4977, 0.153) respectively. (f,h,j,l) Difference between (e,g,i,k) and (b), respectively.

Figure 17 :

 17 Figure 17: Distance evolution during iterations of local search for the inputs in Figure 16 (a) and (b), from the initial transformation in Figure 16 (c), with respect to different depths k = 1, 3, 5, 10 and 15.

Figure 19

 19 Figure19: Distance evolution at each scale level during the multi-scale local search for the inputs in Figure14

  (a) and (b) with the seed (0.365, -0.045, 0.1423) 550 for k = 1, 3, 5, 10, 15.Figure 20 compares between the k-neighbours local search and the multi-scale local search approaches in terms of optimization result of local minimum and the run-time.

Figure 20 :

 20 Figure 20: Comparison between the k-neighbours and the multi-scale local search. Left: difference of distances from the initial DRT to that of local minimum for different k values. Right: total run-time of both approaches for different k values.

Figure 22 :

 22 Figure 22: Results of the local search obtained for the grain numbered 0 in Figure 21, with respect to different k values.

Figure 23 :

 23 Figure 23: Results of the local search obtained for the grains numbered from 1 to 74 in Figure 21, with k = 100.

Table 1 :

 1 Distance results obtained with our method, using the initial transformation from the correlation method [24].

Note that when applying a transformation T v on I s , the transformed digital image I s • T v no longer fits in the N × N support of I s . We deal with this issue by surrounding the image with a padded border. The border values are that of the background. We include this in the distance computation d(I s • T v , I t ).

The algorithm for finding the feasible region R v , and in particular θ min (v) and θ max (v), given 4N 2 tipping surfaces S 1 and S 2 has a complexity O(N 6 )[START_REF] Ngo | On 2D constrained discrete rigid transformations[END_REF]. However, in practice, the initial solution for the gradient descent algorithm is given as a rigid transformation rather than a set of inequalities S 1 and S 2 . Then, the θ min and θ max of R v , as well as its neighbours, can be calculated during the sweeping algorithm. In consequence, knowing θ min and θ max , a representative T v of R v is calculated in constant time using Equation[START_REF] Noblet | 3-D deformable image registration: A topology preservation scheme based on hierarchical deformation models and interval analysis optimization[END_REF].

Note that θ v can be easily computed using Equation[START_REF] Noblet | 3-D deformable image registration: A topology preservation scheme based on hierarchical deformation models and interval analysis optimization[END_REF] with θ min and θ max calculated during the sweeping algorithm.

Let f + ⋆ (resp. f - ⋆ ) be the lowest (resp. uppermost) tipping surface in S + ⋆ (I v ) (resp. S - ⋆ (I v )), that intersects the initial plane at θ v .

6Let S ⋆ be an ordered set containing f + ⋆ (resp. f - ⋆ ) and the next tipping surface that is lower (resp. upper) than f + ⋆ (resp. f - ⋆

). 7

repeat 8 for ⋆ =

1, 2 do 9 Find the next left intersection θ + ⋆ of f + ⋆ (resp. θ - ⋆ of f - ⋆ )with the other surfaces in S ⋆ for θ < θ prev .

repeat

for ⋆ = 1, 2 do

Find the next left intersection θ+ ⋆ of f + ⋆ (resp. θ - ⋆ of f - ⋆ ) with the other

Note that the distance based on signed distance function for grey-level images can be calculated similarly. More precisely, a grey-level image can be modeled by a finite set of its binary level set images. Then, the distance is calculated by summing all the signed distance transforms of the binary level set images.

Acknowledgement

The authors thank to Edward Andò and Gioacchino Viggiani of Laboratoire 3SR, Grenoble, France for providing the data for the experiments, as well as helpful discussions concerning this work. They would like also to thank to Frank Schmidt of University of Freiburg, Germany for helpful discussions and providing the reference [9].

✩ The research leading to these results has received funding from the French Agence Nationale de la Recherche (Grant Agreement ANR-2010-BLAN-0205) and the Programme d'Investissements d'Avenir (LabEx Bézout, ANR-10-LABX-58).

surfaces in S ⋆ for θ < θ prev .

Exchange their order in S ⋆ .

else