É Grossiord

H Talbot

N Passat

M Meignan

P Tervé

L Najman

HIERARCHIES AND SHAPE-SPACE FOR PET IMAGE SEGMENTATION

Keywords: Positron Emission Tomography, segmentation, component-tree, shaping, mathematical morphology, lymphoma

Positron Emission Tomography (PET) image segmentation is essential for detecting lesions and quantifying their metabolic activity. Due to the spatial and spectral properties of PET images, most methods rely on intensity-based strategies. Recent methods also propose to integrate anatomical priors to improve the segmentation process. In this article, we show how the hierarchical approaches proposed in mathematical morphology can efficiently handle these different strategies. Our contribution is twofold. First, we present the component-tree as a relevant data-structure for developing interactive, real-time, intensity-based segmentation of PET images. Second, we prove that thanks to the recent concept of shaping, we can efficiently involve a priori knowledge for lesion segmentation, while preserving the good properties of component-tree segmentation. Preliminary experiments on synthetic and real PET images of lymphoma demonstrate the relevance of our approach.

INTRODUCTION

Positron Emission Tomography (PET) using 18-fluorine fluorodeoxyglucose ( 18 F-FDG) is recognized as the modality of choice for lymphoma imaging, due to its high sensitivity and specificity. PET imaging is now routinely used, not only to detect tumor lesions, but also to assess their metabolic activity, allowing diagnosis, staging and treatment response evaluation. Thus, efficient PET image analysis is essential in clinical practice. However, tumor segmentation is challenging, considering the limitations the modality suffers from: poor spatial resolution compared to anatomical imaging; partial volume effects and intrinsic noise that lowers the quality of reconstructed images; and physiological artifacts. Besides, the segmentation of lymphoma is a difficult task because of multiple tumors location, intensity distribution and contrast to its surrounding tissues. As a consequence, the development of efficient and easy-to-use image processing and analysis tools for guiding the expert's diagnosis is highly relevant.

In contrast to other standard 3D imaging modalities, PET provides metabolic activity visualization, characterized by the intensity of an injected radiotracer. Based on this assertion, most approaches proposed to analyse PET images only rely on this intensity information. On the one hand for metabolic assessment, intensity-based criteria such as the Standardized Uptake Value (SUV) [START_REF] Huang | Anatomy of SUV[END_REF] have been considered. On the other hand, for the detection and segmentation of lesions, intensity-based approaches have been developed. Indeed, the segmentation methods designed during the last ten years were mostly based on thresholding [START_REF] Nestle | Comparison of different methods for delineation of 18 F-FDG PETpositive tissue for target volume definition in radiotherapy of patients with non-small cell lung cancer[END_REF], region-growing [START_REF] Li | A novel PET tumor delineation method based on adaptive region-growing and dual-front active contours[END_REF], classification (FCM [START_REF] Zaidi | Fuzzy clustering-based segmented attenuation correction in wholebody PET imaging[END_REF], FLAB [START_REF] Hatt | A fuzzy locally adaptive Bayesian segmentation approach for volume determination in PET[END_REF]), watershed, or basic mathematical morphology pipelines (see [START_REF] Shepherd | Comparative study with new accuracy metrics for target volume contouring in PET image guided radiation therapy[END_REF] for a recent survey). Practically, such methods generally lead to interactive tools in clinical routines, where the expert user provides regions of interest (e.g., bounding-boxes) and/or seeds, and tune threshold values in order to delineate lesions in the chosen area(s). Despite the low quality of PET images, and based on the assertion that the intensity provides indirect spatial information on organs or lesions gathering the sought metabolic activity, few strategies were designed to integrate anatomical priors in order to improve the segmentation process [START_REF] Wojak | Joint variational segmentation of CT-PET data for tumoral lesions[END_REF][START_REF] Bagci | Joint segmentation of anatomical and functional images: Applications in quantification of lesions from PET, PET-CT, MRI-PET, and MRI-PET-CT images[END_REF]. In contrast to the above intensity-based methods, they generally consider multimodal images (e.g., PET/CT or MRI/PET) to collect supplementary anatomical information, and rely on non-interactive and time-consuming strategies (e.g., deformable models, optimization schemes).

Our purpose is to propose a segmentation strategy combining the advantages of both approaches: the ability to rely on simple and real-time interaction paradigms, based on intensity and thresholding; and the ability to embed anatomical priors that improve the accuracy of segmentation, in particular by discriminating actual lesions from organs fixating the radiotracer. To reach that goal, we propose for the first time to consider the morphological hierarchy paradigm [START_REF]Mathematical Morphology: From Theory to Applications[END_REF]Ch. 7] for designing PET-oriented connected operators. In Sec. 2, we propose to use the notion of component-tree [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF] as a hierarchical data-structure allowing us to perform PET segmentation. Indeed, the component-tree is well-fitted, both for handling 3D images where the structures of interest correspond to extremal intensity values (as already proved in angiographic imaging [START_REF] Dufour | Filtering and segmentation of 3D angiographic data: Advances based on mathematical morphology[END_REF]) and for developing interactive and real-time segmentation methods [START_REF] Passat | Interactive segmentation based on component-trees[END_REF]. In Sec. 3, we describe a strategy based on the recently introduced notion of the shapings [START_REF] Xu | Morphological filtering in shape spaces: Applications using tree-based image representations[END_REF], that consists of using a two-layer component-tree, namely a "tree of tree", to embed additional information such as geometrical priors. Since the induced data-structure remains a tree, the desired properties described in Sec. 2, namely user interaction, thresholding paradigm, and real-time computation, remain valid, while naturally handling higher-level knowledge. In Sec. 4, we validate this methology on phantom images and on 3D PET images of lymphoma, where results are compared with expert's manual delineation.

COMPONENT-TREES AND GENERALIZED THRESHOLDING

Component-trees

From a structural point of view, the component-tree [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF] of an image I is a data-structure T that models all the connected components (i.e., the maximally connected regions) of I, obtained from its successive binary level (thresholded) sets. More formally, the component-tree T is composed of a set of nodes Θ, namely the connected components of the binary level sets, and of edges between these nodes, that correspond to the transitive opening of the partial order relation of inclusion ⊆ on Θ. Each node K ∈ Θ is associated to the (maximal) value v of the level set where the associated connected component appears.

Various algorithms exist for computing the component-tree of an image, and the most efficient run in quasi-linear time [START_REF] Carlinet | A comparative review of component tree computation algorithms[END_REF]. In addition, retrieving an image from its component-tree is a trivial process that can be carried out in linear time. Processing an image via its component-tree is then a low-cost operation.

Moreover, the component-tree is a lossless image model. Indeed, if we consider the image I in its functional form, i.e., as a mapping I : Ω → V , where Ω is the image support, i.e., the set of its voxels, and V is its set of grey-level values, then the image I can be expressed as the supremum of the nodes of T :

I = K∈Θ CK (1)
where CK : Ω → V is the cylinder function such that CK (x) = v for any x ∈ K with v the maximal value that defines K as a connected component of the level set of I at threshold v, while CK (x) = 0 anywhere else. The composition formula (1) leads to a well-defined image for Θ, but also for any subset of nodes of Θ. In other words, it is possible to filter the image I by discarding some of the nodes, and then reconstructing a resulting image from the preserved nodes. Each point x ∈ Ω in the filtered image then presents a value that is lower or equal to the initial image; the induced operators are then antiextensive.

This anti-extensive filtering paradigm was formalized in [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF][START_REF] Jones | Connected filtering and segmentation using component trees[END_REF], and is basically composed of three successive steps:

(i) construction of the component-tree T associated to I;

(ii) reduction of Θ, leading to a reduced component-tree T ; and

(iii) reconstruction of a filtered image I induced from T .

I (i) -----→ T     (ii) I ≤ I (iii) ← -----T (2)

PET image segmentation with component-trees

In the case of PET, the component-tree of an image is organized from the lowest to the highest values. More precisely, for a given PET image I, the root of its component-tree T , i.e., the node K that is equal to the whole support Ω, is obtained from the lowest level set, at value 0. Conversely, the leaves of T , i.e., the nodes K that are minimal elements for the inclusion relation ⊆, are the flat zones of the image of (locally) maximal values, and then correspond to high metabolic activity areas of lesions and hyperfixating organs.

The main strength of the component-tree, is that it intrinsically models the space of all the potential thresholding operations that can be carried out on an image. In addition, the choice of a solution within this data-structure can be performed in linear time, by scanning the (decreasing) grey-level values along the branches of the tree, from its leaves to its root. In particular, this process is compliant with all major strategies, namely fixed thresholding (horizontal cut in the tree), adaptive thresholding (branch-specific cut in the tree), or bounding-box thresholding (explicit choice of leaves or branches to cut in the tree).

Component-trees then constitute a highly relevant model for intensity-based analysis of PET images. Its low space and time cost properties are crucial to perform PET image segmentation. In the next section, we propose a framework for performing componenttree-based segmentation by embedding a priori geometric information, without losing these complexity properties.

FROM INTENSITY TO ATTRIBUTE SEGMENTATION:

THE SHAPING PARADIGM

Shaping: Trees of trees

Another strength of component-trees, that derives from the decomposition of the image into elementary patterns both in terms of space and intensity, is the ability to associate an attribute a ∈ A to each node of the tree, that models additional information. This attribute-based approach has been considered in the antiextensive filtering framework described in Diag. [START_REF] Nestle | Comparison of different methods for delineation of 18 F-FDG PETpositive tissue for target volume definition in radiotherapy of patients with non-small cell lung cancer[END_REF]. The choice of the nodes to preserve or discard depends on a Boolean predicate ρ : A → {true, false} acting on the attributes. When the attribute is increasing (i.e., ρ(K) is true for a given node K implies the same for all the nodes between K and the root), the attribute can be handled similarly to the grey-level of the nodes, in the sense where the node selection proceeds as a simple thresholding. For instance, this is the case for attributes such as the number of voxels of the node, used for area/volume denoising filters. However, most attributes are not increasing, in general. This means that the validity of the predicate for a given node does not imply its validity for the rest of the branch. This is not an issue in the context of automated processing, where the tree is fully scanned, and ad hoc pruning policies can be applied a posteriori to homogeneize the results in the different branches of the tree (see [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF][START_REF] Jones | Connected filtering and segmentation using component trees[END_REF] for a discussion on this topic).

In our case, as we wish to carry out threshold-based interactive segmentation, working with increasing attributes is crucial. To tackle this issue, we propose to rely on the recently introduced notion of shapings [START_REF] Xu | Morphological filtering in shape spaces: Applications using tree-based image representations[END_REF]. Basically, shapings consists of considering the component-tree T as an image. Indeed, from a functional point of view, it can be defined as a mapping T : Θ → A, where voxels are replaced by nodes, while grey-levels are replaced by values of the attribute. This approach is tractable only if the space of the attribute values A is equipped with a total order relation, i.e., can be modeled as (a subset of) R or Z; this is however the case of most attributes, in particular numerical ones.

In this case, the "image" T can be thresholded, and for any level set, the connected components, i.e., the sets of nodes that are linked by edges in T , can be computed and ordered with respect to the inclusion relation. It is then possible to build a component-tree of the component-tree T . This "tree of tree", noted T T , can be processed as any other component-tree, and we can perform anti-extensive filtering with it. The virtue of this new tree is that the attribute computed from the nodes of T is now increasing in T T , which can then be processed by a thresholding approach. The overall procedure remains quasi-linear in time and space, since we only duplicate the standard component-tree anti-extensive filtering process, as illustrated in Diag. (3):

I (i) -----→ T (i) -----→ T T   (ii)   I ≤ I (iii) ← -----T (iii) ← -----T T (3)

Attribute-based segmentation of PET images

We consider the paradigm of shapings to develop an interactive attribute-based procedure for discriminating active lesions from hyperfixating organs in PET images, in the context of lymphoma. As stated above, an intensity-based strategy, relying either on a single component-tree or any other -classification, region-growing, watershed -intensity-based approach, may over-emphasise the high intensity areas, including true positives, but also many false positives.

Our working hypothesis, that is relevant in the case of lymphoma -where the number of lesions is high and their shape is compact, compared to other fixating organs (e.g., heart, kidneys, etc.) -is that a compacity/elongation criterion can allow us to better discriminate the active lesions.

Then, we propose to define an attribute E : 2 Ω → R that associates to any connected subset of Ω, a numerical value describing how the shape is compact, i.e., similar to a sphere, or elongated. This attribute is defined by computing the eigenvalues λ1 ≥ λ2 ≥ λ3 associated to the eigenvectors of the 3 × 3 matrix of inertia of the (binary) shapes K ⊆ Ω. We define

E(K) = λ3/λ1 (4)
This is the ratio between the lowest and highest eigenvalues, lying in [0, 1]. If E(K) is close to 1, then λ3 and λ1 (and thus λ2) present similar values, and K then has a compact shape. When E(K) decreases toward 0, K becomes more elongated. The PET image is first processed to build its component-tree T , based on the total order of the metabolic signal. As stated above, in this first tree, the nodes that are close to the leaves correspond to the areas of highest intensity, but can be either lesions or hyperfixating organs. For each node K ∈ Θ of this tree, the attribute value E(K) ∈ [0, 1] is then computed and stored at the node.

The second component-tree T T is then defined from the "image" T : Θ → [0, 1], that takes its values in the shape space induced by E. In this new tree, the root is the node of maximal value, while the leaves are those of locally minimal values. This second tree T T is a min-tree (i.e., a component-tree based on (R, ≥)), while T was a max-tree (i.e., a component-tree based on (R, ≤)). The nodes that are close to the leaves then correspond to the more elongated areas, i.e., the sets of connected components of Ω that are of least relevance, and should then be removed preferentially.

In order to segment the active lesions versus the hyperfixating organs, we propose two distinct strategies. The first is a thresholding along the branches of the tree T T , thus pruning the distal parts of the tree, composed of the nodes that have a compactness value E lower than a chosen value C ∈ [0, 1]. The second is a differential thresholding that considers the "gap" between the compactness value EK of the node K and the value EL of the leave of its branch. This gap EK -EL still varies in [0, 1], with 0 values at the leaves. This criterion remains increasing, thus authorising a (relative) thresholding approach similar to the first strategy, but with a finer behaviour. This latter strategy is considered for the experiments described in the next section.

EXPERIMENTS AND RESULTS

Phantom images

To validate our approach, we first used the NEMA 2007 IEC image phantom. This phantom contains six spherical "lesions" with diameters 10, 13, 17, 22, 28 and 37 mm. The spheres were filled with varying 18 F concentrations to three different signal to background ratios (SBRs), namely 2.59:1, 5.06:1 and 8.40:1. Data were acquired on a Philips GEMINI TF Time-of-Flight PET/CT scanner. The PET images were obtained with the 3D LOR-RAMLA reconstruction and CT-based attenuation correction. They present a spatial resolution of 2.0 mm (Fig. 1, first row).

For comparison, all images were registered to the highest contrast image. Reference segmentation was obtained by closest interactive thresholding to the known geometry on this image. The ROC curves of Fig. 2 measure the set differences between the reference segmentation and our shapings result. It illustrates the impact of the signal loss on the segmentation results. The obtained results are near perfect for the image with a very high SBR of 8.40:1. However, even for more realistic SBRs, including that of value 2.59:1, which is low with respect to real cases, the obtained results remain satisfactory, thus underlining the relevance of the shapings approach.

Real images

Finally, we qualitatively evaluated our method on a series of 13 PET images of lymphoma patients presenting significant FDG-PET tumor radioactivity uptake before radiotherapy. The acquisitions were carried out on several PET/CT scanners (Philips ALLEGRO and GEMINI DUAL GS, Siemens Biograph Sensation 16, GE DISCOV-ERY ST) with their associated reconstruction algorithms (RAMLA, TF MLEM, OSEM). They were acquired following standard protocol for PET cancer imaging: PET acquisition using a scan time of 3 minutes per bed position, between the base of the skull and mid-tigh, one hour after the peripheral intravenous injection of an 18 F-FDG dose of 4 to 6 MBq/kg in patients, fasted for at least 6 hours before FDG injection. PET image were obtained with CT-based attenuation correction and their spatial resolution (FWHM) varied from 4×4×2 mm 3 to 4 × 4 × 4 mm 3 .

Our approach was applied to these images and compared to global thresholding. An example of result is illustrated in Fig. 3. We observe that the proposed method improves the robustness to signal heterogeneity in the image; this fact derives from the component-tree based approach, that intrinsically provides adaptive thresholding. The shapings paradigm also allows us to take into account geometric priors, improving in particular the robustness to texture effects, for instance here in the right lung. Furthermore, shapings share with the global thresholding the advantages of low parametricity, easy interaction and real-time computing.

CONCLUSION AND PERSPECTIVES

We have introduced the first application of morphological hierarchies for the segmentation of PET images. We have shown that the component-tree handles the principal approaches of intensity-based segmentation, subsuming the thresholding strategies considered until now. We have also proposed a solution based on the shapings paradigm to embed anatomical priors, while preserving a thresholdbased strategy, that allows us to satisfy important properties such as real-time computing, and user-friendly interaction. For future work, we plan to use several attributes, to enrich the anatomical priors guiding the segmentation. Second, we also plan to use morphological (CT, MRI) data, in addition to PET images, as already proposed in recent works [START_REF] Wojak | Joint variational segmentation of CT-PET data for tumoral lesions[END_REF][START_REF] Bagci | Joint segmentation of anatomical and functional images: Applications in quantification of lesions from PET, PET-CT, MRI-PET, and MRI-PET-CT images[END_REF]. These are indeed prerequisites to retrieve accurate anatomical information, both from geometrical and relational points of view. As a side effect, we will need to handle vectorial attributes (resp. multimodal imaging), and thus to define the attribute space of T (resp. I) in a partial order framework, which are not handled by component-trees. To tackle this issue, we will cross the con- cept of shapings with the recently introduced notion of componentgraph [START_REF] Passat | Component-trees and multivalued images: Structural properties[END_REF], in order to develop the next generation of hierarchical approaches for PET image segmentation.

1 Fig. 1 .

 11 Fig. 1. First row: phantom PET images, with three different SBRs. Second row: segmentation results. MIP visualization. See Sec. 4.1.

Fig. 2 .

 2 Fig.2. ROC curves for the relative thresholding based on attributebased approach, for various SBRs (see Fig.1).

Fig. 3 .

 3 Fig. 3. Left: PET image. Middle: optimal global thresholding. Right: relative shaping segmentation. Maximum intensity projection. See Sec. 4.2.

ACKNOWLEDGEMENTS

The phantoms used in section 4.1 were provided courtesy of Irène Buvat from Université Paris 7 & Université Paris 11 and Marine Soret, Hôpital Inter Armées du Val de Grâce, Paris.


Phuc Ngo

Yukiko Kenmochi

Nicolas Passat

Hugues Talbot

Sufficient conditions for topological invariance of 2D images under rigid transformations ⋆

Keywords: 2D digital image, discrete rigid transformation, topology, simple point, DRT graph, Eulerian model

In R 2 , rigid transformations are topology-preserving operations. However, this property is generally no longer true when considering digital images instead of continuous ones, due to digitization effects. In this article, we investigate this issue by studying discrete rigid transformations (DRTs) on Z 2 . More precisely, we define conditions under which digital images preserve their topological properties under any arbitrary DRTs. Based on the recently introduced notion of DRT graph and the classical notion of simple point, we first identify a family of local patterns that authorize topological invariance under DRTs. These patterns are then involved in a local analysis process that guarantees topological invariance of whole digital images in linear time.

Introduction

In 2D, rigid transformations (i.e., rotations composed with translations) are involved in numerous image processing/analysis tasks, e.g., registration [START_REF] Zitová | Image registration methods: A survey[END_REF] or tracking [START_REF] Yilmaz | Object tracking: A survey[END_REF]. In such applications, the images are generally digital, and can then be considered as functions I : S → F from a finite subset S ⊂ Z 2 to a value space F. While rigid transformations are topology-preserving operations in R 2 , this property is generally lost in Z 2 , due to the discontinuities induced by the mandatory digitization from R to Z. In particular, discrete rigid transformations (DRTs) -that include discrete rotations [START_REF] Jacob | On discrete rotations[END_REF][START_REF] Andres | The quasi-shear rotation[END_REF][START_REF] Nouvel | Configurations induced by discrete rotations: Periodicity and quasiperiodicity properties[END_REF][START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]-are not guaranteed to preserve the homotopy type of digital images, as exemplified in Fig. 1.

In this article, we study this specific issue. More precisely, we investigate some conditions under which digital images preserve their topological properties under any arbitrary DRTs, by considering the Eulerian (i.e., backwards) transformation model. To reach this goal, we consider (i) the notion of DRT graph, recently introduced by the authors in [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF][START_REF] Ngo | Combinatorial properties of 2D discrete rigid transformations under pixel-invariance constraints[END_REF], that defines a combinatorial model of all the rigid transformations of a digital image, and (ii) the classical notion of simple point [START_REF] Kong | Digital topology: Introduction and survey[END_REF][START_REF] Couprie | New characterizations of simple points in 2D, 3D, and 4D discrete spaces[END_REF], that provides sufficient conditions to guarantee the preservation of homotopy type.

By combining these two notions, we first propose a way to determine transformed images which have the same homotopy type as the initial one, by scanning the whole DRT graph associated to this image. Then, we show that this global approach, which presents a polynomial complexity, can be simplified into a local approach, based on a Fig. 1. Left: a binary digital image and the grid modeling its discrete structure. Middle: a rigid transformation applied on this grid. Right: the resulting transformed image, with a homotopy different from the initial one (the black pixels, in the 8-adjacency, have been split).

spatial decomposition of the image into covering samples. In order to do so, we identify a family of local patterns that authorize topological invariance under DRTs. These patterns can then be involved in a procedure based on look-up tables (LUT) that guarantee topological invariance of a whole digital image in linear time.

The article is organised as follows. Sec. 2 presents background notions related to rigid transformations and digital topology. Sec. 3 describes the topological issues induced by DRTs. Sec. 4 explains how DRT graphs and simple points can be combined to evaluate topological invariance under DRTs, leading to an algorithm detailed in Sec. 5. Experiments are proposed in Sec. 6, while Sec. 7 concludes the article.

Background notions

(Discrete) rigid transformations

In R 2 , a rigid transformation (i.e., a transformation composed of a translation and a rotation) is expressed as a bijection T : R 2 → R 2 defined, for any x = (x, y) ∈ R 2 by

T (x) = cos θ -sin θ sin θ cos θ x y + a b with a, b ∈ R and θ ∈ [0, 2π[ (1)
Such a transformation (also noted T abθ ) is unambiguously modeled by the triplet of parameters (a, b, θ). It is not possible to apply directly T on a digital image I : S → F, since there is no guarantee that T (x) ∈ Z 2 , for any x ∈ S ⊂ Z 2 . The handling of discrete rigid transformations (DRTs) then requires the definition of a function T : Z 2 → Z 2 , which is the "discrete analogue" of T . Considering the standard rounding function D : R 2 → Z 2 , this can be conveniently performed by setting T = D • T , as illustrated on the diagram below.

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

T ( p) = D • T ( p) = [p cos θ -q sin θ + a] [p sin θ + q cos θ + b] (2)
In general, this function is not bijective. However, by setting

T -1 = D • T -1 : Z 2 → Z 2 ,
we can define the transformed digital image

I • T -1 : Z 2 → F with respect to T . Note that T -1 is not the inverse function of T in general. t z y x Fig. 2.
Examples of simple points (x, y) and non-simple points (z, t). Modifying the value of z would merge two black connected components, while modifying the value of t would create a white connected component. In both cases, the homotopy type of the image would be modified.

Digital topology

Several frameworks are available to model the topological structure of a digital image. In Z 2 , most of these frameworks (see, e.g., [START_REF] Khalimsky | Topological structures in computer science[END_REF][START_REF] Kovalevsky | Finite topology as applied to image analysis[END_REF]) can be conveniently unified within the frequently used -and also simple-framework of digital topology [START_REF] Kong | Digital topology: Introduction and survey[END_REF]. In this framework, the topological notions derive from a graph structure induced by two adjacency relations, namely the 4-and 8-adjacencies, which are defined for any two points p, q ∈ Z 2 such that p and q are 4-adjacent (resp. 8-adjacent

) if || p -q|| 1 ≤ 1 (resp. if || p -q|| ∞ ≤ 1)
. It is well known that, to deal with topological paradoxes related to the digital version of the Jordan theorem, we generally use in one binary digital image a pair of different adjacency relations, and denote as (α, β) where α and β are adjacency relations for foreground (black) and background (white) pixels respectively. In 2D, we consider in particular (α, β) = (4,8) or [START_REF] Ngo | Combinatorial properties of 2D discrete rigid transformations under pixel-invariance constraints[END_REF][START_REF] Andres | The quasi-shear rotation[END_REF].

In the graph-based framework of digital topology, the concept of simple point [9,10] (see Fig. 2) relies on the local notion of adjacency and on the induced global notion of connectedness. The simple points provide a way to characterise the preservation of topological properties in a (binary) image during its transformation. Practically, a pixel x ∈ S of an image I : S → F is simple if its binary value can be switched without modifying the topological properties of I. In particular, the simplicity of a pixel can be tested, in constant time, by only studying its 3 × 3 neighbourhood [START_REF] Kong | Digital topology: Introduction and survey[END_REF]. We will say that two images I and I ′ are simple-equivalent [START_REF] Bertrand | A note on 3-D simple points and simple-equivalence[END_REF] if I ′ is obtained from I by iteratively modifying (successive) simple points. Thus I and I ′ present the same homotopy type.

Discrete rigid transformations: topological issues

A (continuous) rigid transformation T establishes a bijection from R 2 to itself. By opposition, due to the digitisation process induced by D (see Eq. ( 2)), a discrete rigid transformation T is, most of the time, not a bijection from Z 2 to itself.

It is plain that for any three distinct pixels

x 1 , x 2 , x 3 ∈ Z 2 , we have max i, j∈{1,2,3} {||x i - x j || 2 } ≥ √ 2.
This leads to the following definition that enables to characterise the status of a pixel; there are only three possibilities, as illustrated in Fig. 3.

Definition 1 For a pixel x ∈ Z 2 and a given DRT T , let M(x) = {y ∈ Z 2 | T (y) = x}.
-If |M(x)| = 0, we say that x is a null pixel.

-If |M(x)| = 1, we say that x is a single pixel.

-If |M(x)| = 2, we say that x is a double pixel. Similar notions for the case of discrete rotations can be found in [START_REF] Nouvel | Configurations induced by discrete rotations: Periodicity and quasiperiodicity properties[END_REF][START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF].

In particular, a discrete rigid transformation T behaves like a bijection for single pixels. However the possible existence of null (resp. double) pixels may forbid T to be a surjection (resp. an injection). Null and double pixels thus raise topological issues in both Lagrangian and Eulerian transformation models (see Sec. 3.1). In addition to these "cardinality-based" issues, supplementary topological problems are induced by the alteration of adjacency relations between pixels (see Sec. 3.2).

Transformation models

Two standard transformation models can be considered: the Lagrangian (or forwards) and the Eulerian (or backwards) models. The Lagrangian model consists of observing T (x) for every pixel x in the initial space, while the Eulerian model consists of observing T -1 (y) for every pixel y in the transformed space. These models are equivalent in R 2 , since T is bijective. In Z 2 , they are however generally distinct, since T and T -1 may not be inverse functions.

Depending on each model, null and double pixels lead to different interpretations. In the Lagrangian model (see Fig. 4(a)), a double pixel in the transformed space may receive two different pixel values, and a null pixel receive no pixel value. While this may be conveniently handled in the case of binary images (which can be considered in a set-based paradigm), it can lead to correctness and completeness issues in the case of multivalued images. In the Eulerian model (see Fig. 4(b)), a double pixel of the initial space may transfer its value to two pixels of the transformed space, while the value of a null pixel will be lost.

In this preliminary work, we consider the Eulerian model which enables us to focus on the topological issues raised by the alteration of adjacency relations (see Sec. 3.2), and as the Lagrangian model is fraught with additional difficulties (see Sec. 7). For the sake of readability, our study is carried out in the context of binary images, but the introduced methodology remains valid in the case of multivalued images (see Sec. 7).

Adjacency alterations

In order to illustrate the topological issues raised by the alterations of adjacency relations during discrete rigid transformations, let us consider a 2 × 2 pixel sample of the transformed space (see Fig. 5 sample is generated (see Fig. 5(b)). Despite local adjacency alterations between pixels, the global topology of the sample may sometimes be preserved (see Fig. 5(c)). Unfortunately, such local alterations may also lead to topological alterations in the sample (see Fig. 5(d)), and further in the whole image possibly. In the next section, we propose an algorithm enabling the detection of potential topological changes during a DRT. On the contrary, this algorithm can be used to guarantee the topological invariance between an image and all of its transformed ones. This algorithm is based on (i) the recently introduced notion of DRT graph [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF][START_REF] Ngo | Combinatorial properties of 2D discrete rigid transformations under pixel-invariance constraints[END_REF], and (ii) the classical notion of simple point [START_REF] Kong | Digital topology: Introduction and survey[END_REF][START_REF] Couprie | New characterizations of simple points in 2D, 3D, and 4D discrete spaces[END_REF]. The first notion provides a way to exhaustively explore the space of transformed (sub)images while the second provides information on the possible topological modifications when performing such an exploration.

Mathematical tools for topological verification of images under rigid transformations 4.1 Discrete rigid transformation graph

In opposition to rigid transformations in R 2 , discrete rigid transformations (DRTs) are not continuously defined w.r.t. the parameters a, b and θ. In particular, the parameter 2)). In fact, each 2D surface is corresponded to an elementary modification of one pixel.

From a theoretical point of view, each 3D open cell can be seen as the equivalent class of the rigid transformations T that leads to a same DRT T = T • D whose boundaries are the 2D surfaces. By mapping any 3D cell onto a 0D point and any 2D surface onto a 1D edge, the combinatorial structure of the parameter space can be modeled in a dual way by a connected graph, as illustrated in Fig. 6. In particular, each 0D point corresponds to a transformed image generated by the associated T and an 1D edge between two cells indicates that the two associated images differ at exactly one pixel.

Definition 2 (DRT graph [7]) A DRT graph G = (V, E) is defined such that:

-any vertex v ∈ V models a 3D open cell and associates to a transformed image; -any (labeled) edge e = (v, w, ( p, p ′ )) ∈ E models a 2D discontinuity surface between the transformed images corresponding to the DRTs v and w which differ at the single pixel p ′ . Note that p is the pixel corresponding to p ′ in the original image.

The label ( p, p ′ ) on each edge is -implicitly-associated to a function indicating the value modification of the pixel p ′ that differs between the DRTs v and w. More precisely, the value of p at the vertex v is defined by I v ( p ′ ) = I( p) where I : S → F is the original image function. After the elementary rigid motion at e, we then obtain a new image I w by simply changing the pixel value at p as I w ( p ′ ) = I( p + δ) where δ = (±1, 0) or (0, ±1). Note that δ corresponds to an elementary motion, i.e., a smallest pixel movement, that changes either x-or y-coordinate by 1.

It was proved in [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF] that the DRT graph associated to a digital image of size N 2 has a space complexity of O(N 9 ) (and can be built with a similar time complexity [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF][START_REF] Ngo | Combinatorial properties of 2D discrete rigid transformations under pixel-invariance constraints[END_REF]). Note that the structure of the DRT graph depends only on the support of the given images, but not on their pixel values. By construction, the DRT graph provides all the transformed images of a given image I. In particular, these transformed images can be generated by progressively and exhaustively scanning the DRT graph.

Remark 3 Let I : S → F, and G = (V, E) be its DRT graph. For each edge e = (v, w, ( p, p ′ )) (i.e., each elementary modification of a pixel p ′ ∈ S ), two cases can occur:

(i) I v ( p ′ ) = I w ( p ′ ), i.e., the images of I by the DRTs v and w are equal (I v = I w );

(ii) I v ( p ′ ) I w ( p ′ ), i.e., I v I w . In the (considered) case of binary images, the value of p ′ may then be flipped from white to black (or vice versa), and this may constitute the only modification between the transformed images of I by the DRTs v and w.

Topological analysis of binary images from DRT graphs

From a DRT graph, one can generate exhaustively all the DRT images of an image I. Moreover, from Rem. 3, we know that this can be done by modifying (at most) one pixel value between two successive transformed images. In the case of binary images, it is actually possible to check if such a local modification involves a simple point.

Practically, the edges of a DRT graph G = (V, E) can then be classified in two categories: those that do not modify the topology of the transformed images (i.e., the edges that correspond to case (i) in Rem. 3, as well as those that correspond to case (ii) for which p is a simple point); and those that modify this topology (i.e., the edges that correspond to case (ii) in Rem. 3 for which p is not simple).

The partial graph G ′ = (V, E ′ ) is obtained by maintaining only the edges E ′ ⊆ E of the first category. G ′ is composed of connected components of vertices whose associated transformed images are simple-equivalent (see Sec. 2.2), and thus have the same homotopy type. In particular, the connected component contains the vertex u corresponding to the initial image I, as well as those corresponding to transformed images obtained from I by elementary motion sequences which are topology-preserving. This specific set of vertices can be straightforwardly computed by using a standard spanningtree algorithm, initialized from the seed vertex u (see Alg. 1).

Remark 4 The connected component of G ′ that contains u may constitute only a strict subset of the vertices/transformed images that are simple-equivalent to u/I. Indeed, the edges of the DRT graph G only model the local modifications associated to DRTs. In particular, there may exist other series of local modifications relying on simple points but not modeled in the DRT graph. In other words, the analysis of the DRT graph provides sufficient (but not necessary) conditions to guarantee homotopy-type preservation.

In the case where V ′ = V (see Alg. 1), i.e., when all the vertices of the DRT can be reached from u by a sequence of edges involving only simple points, the algorithm successfully detects -as a side effect-an image I that is actually topologically invariant under any DRTs. The algorithmic cost of this algorithm is directly linked to the size of the DRT graph, that is O(N 9 ). This algorithmic complexity is indeed reached in the worst cases. In the next section, we show that this problem can however be decomposed spatially, thus leading to a practical, lower complexity algorithm.

Local evaluation of topological invariance under DRTs

In the previous section, we have proposed to explore the whole DRT graph of an image I in order to check its potential topological invariance under DRTs. For each edge e = Algorithm 1: Construction of simple-equivalent DRT images ,w,( p,p ′ )) of the DRT graph, we verify that the pixel p ′ whose value is modified at this edge is actually a simple point for the transformed images I v and I w . This test is performed locally, more precisely in the 3 × 3 neighbourhood centered on p ′ in the transformed space.

Input: A DRT graph G = (V, E); the vertex u ∈ V associated to I. Output: A connected partial subgraph G ′′ = (V ′ , E ′′ ) of vertices simple-equivalent to u/I. 1 V ′ ← {u}; E ′′ ← ∅; S ← {u}; U ← V; 2 while S ∅ do 3 Let v ∈ S; S ← S \ {v}; 4 if v ∈ U then 5 U ← U \ {v}; 6 foreach e = (v, w, (p, p ′ )) ∈ E such that w ∈ U do 7 if p ′ is a simple point then V ′ ← V ′ ∪ {w}; E ′′ ← E ′′ ∪ {e}; S ← S ∪ {w}; (v
We now take advantage of the local nature of these tests to develop a spatial decomposition strategy that will lead to a local version of the previously proposed global method. To this end, we first need to introduce basic notions and properties related to the influence of DRTs on pixel neighbourhoods.

Neighbourhoods and DRTs

Let p ∈ Z 2 be a pixel. We define the neighbourhoods of p as follows:

N 8 ( p) = {q ∈ Z 2 | ||q -p|| 2 < 2}; N 20 ( p) = {q ∈ Z 2 | ||q -p|| 2 < 2 √ 2}.
The first 3 × 3 neighbourhood is classically used in digital geometry and topology. The second corresponds to a 5 × 5 square from which the 4 extremal corner pixels have been removed. We provide the following property where we consider any arbitrary DRT T :

Z 2 → Z 2 .
Property 5 Let p ∈ Z 2 and q ∈ N 8 ( p). We then have T -1 (q) ∈ N 20 (T -1 ( p)).

Proof This property derives from the above definitions of the neighbourhoods, and from the fact that a DRT T (due to the digitization induced by D, see Eq. ( 2)) implies a possible (strict) majoration of √ 2 for the distance between transformed points, w.r.t. its associated rigid transformation T .

A local approach for topological analysis

As stated above, a DRT graph models all the rigid transformations of a given digital image I. Despite the fact that the space of these transformations is actually infinite, the DRT graph is defined as a finite structure. This restriction can be made without loss of correctness/completeness by considering translation invariance. Indeed, a rigid transformation is defined as a composition of a rotation and a translation (Eq. ( 1)). In particular, a rigid transformation T = t • r, composed of a rotation r around the origin and a translation t ∈ R 2 , can be also expressed as

T = t • t ′-1 • r ′ • t ′ ,
where t ′ is the translation by a vector p ∈ Z 2 , and r ′ is the rotation around p. Such a translation in Z 2 (that induces no topological modification, since the whole image is translated), allows any pixel p of S of the image I to be considered as the origin of S .

Let come back to the DRT graph G = (V, E) considered until now, and involved in the global process defined in the previous section. We now focus on an edge e = (v, w, ( p, p ′ )) of G. Obviously, the vertices v and w exist in G, and also in the (equivalent) DRT graph where we consider p as origin. In the later graph, any edge that does not involve in its label a pixel of N 20 ( p) has no influence on the topological modifications in N 8 ( p ′ ) (see Prop. 5). Without loss of correctness, such topological modifications at p ′ in N 8 ( p ′ ) (and thus of p ′ in the whole image) only depends on the part of the DRT graph that corresponds to the restriction of I to N 20 ( p), denoted by I |N 20 (p) . Based on these considerations, it appears that if for any pixel q in the initial image I, the restriction I |N 20 (q) does not lead to topological modifications under any DRTs, then the same conclusion holds for the whole image I. In other words, every elementary topological change occurring on the DRT graph G of I can be observed locally. Therefore, we need only to verify the topological invariance for every pixel of I in its neighbourhood N 20 in the original binary image.

Proposition 6 Given a binary image I : S → F, for every p ∈ S if I |N 20 (p) is a local binary configuration which is topologically invariant under any rigid transformations, then the image I is topologically invariant.

We assume that the modified pixel p ′ at each elementary rigid modification of the DRT graph corresponds to the origin o 2 in the transformed image, and that the corresponding pixel p is the origin o 1 in the original image of size in its N 20 (o 1 ). Then, we simply need to construct the DRT graph with p ′ ∈ [-1, 1] 2 , p ∈ [-2, 1] 2 (i.e., the DRT graph of edges labels ( p, p ′ ) with this constraint), denoted by G p = (V p , E p ). We use Alg. 1 proposed in Sec. 4.2 to verify in G p the topological equivalence between two adjacent vertices v and w whose edge has the label (o 1 , o 2 ). If every edge has topologically equivalent vertices, then the center point o 1 of such a configuration is topologically invariant under any rigid transformations. This approach, in particular, leads to the following consideration: if we study the topological invariance property for all the binary image configurations of N 20 (o 1 ), we can identify a family of elementary configurations that authorise topological invariance under DRTs.

From Prop. 6, we propose a look-up-table-based algorithm for characterizing the topologically invariance property of any binary image. More precisely, we generate a set P 4 (resp. P 8 ) which contains only topologically invariant configurations in (4, 8)-(resp. (8, 4)-) adjacent relations. Then we use P 4 and P 8 to verify whether the given image is topologically invariant. The method for building P 4 and P 8 is given in Alg. 2. Let C be the set of all binary image configurations of size N 20 , which is used to build P to study the topologically invariance property of each configurations in C. We store in P 4 and P 8 the subset of C containing only the topologically invariant configurations w.r.t. the (4, 8)-and (8, 4)-adjacent relations. Using Alg. 2, we obtain sets of 10 643 and 19 446 topologically invariant configurations in P 4 and P 8 respectively. Fig. 7 shows some elements of P 4 and P 8 .

Experiments

In this section, we illustrate the relevance of our approach by presenting images which have been detected as topology-invariant (see Fig. 8(a-c),(d-f)), or topology-variant (see Fig. 8(g,h), (i,j), (k,l)). Thanks to our LUT-based approach, such a detection can be carried out in linear time w.r.t. the image size.

As mentioned in Rem. 4, we only have a sufficient condition for homotopy-type preservation, so far we do not have a proof for a necessary condition. Nonetheless we have not found any example for which our algorithm fails to characterize its topological invariance. We have proposed an algorithmic process for determining the topological invariance of digital images under any discrete rigid transformations. This work is based on the recently introduced notion of DRT graph [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF][START_REF] Ngo | Combinatorial properties of 2D discrete rigid transformations under pixel-invariance constraints[END_REF], which presents a polynomial complexity that generally forbids its practical application on whole images. Nevertheless, DRT graphs have been successfully involved in a preliminary local analysis that finally led to a low complexity methodology, relying on image spatial decomposition.

Beyond its theoretical aspects, this work may contribute to the better understanding of the relationships that exist between geometry and topology in the framework of digital imaging, where both notions are more strongly linked than in continuous spaces.

This study was carried out in the context of binary images. However, it remains relevant whenever a notion of simple point (or more generally a local characterisation of topology preservation) is available. This is verified, for instance, in the context of n-ary images [START_REF] Mazo | Topology on digital label images[END_REF]. On the other hand, only the Eulerian (backwards) model has been considered in this study. In future work, we will extend these results to the case of the Lagrangian (forwards) model. Note that additional difficulties arise in the Lagrangian model, such as double pixels in the transformed space that may receive two different values, and null pixels that do not have any value. The Lagrangian model thus involves a value decision problem for such pixels.

Fig. 3 .

 3 Fig. 3. Left: a digital image support and the grid modeling its discrete structure. Right: examples of a null pixel (in green), a single pixel (in blue) and a double pixel (in red) with respect to a discrete rigid transformation.

Fig. 4 .Fig. 5 .

 45 Fig. 4. The interpretations of double pixels (left figures) and null pixels (right figures) in the context of discrete rigid transformations for (a) the Lagrangian and (b) the Eulerian models. a b c d (a)

Fig. 6 .

 6 Fig. 6. A part of the parameter space subdivided by four 2D surfaces corresponding to the discontinuities of DRTs (left), and the associated part of the DRT graph (right).

4 and P 8 ,

 8 |C| ≤ 2 20 . From Rem. 3 we have I v I w if I v (o 1 ) I w (o 1 ), where I v (o 1 ) = I(o 2 ) and I w (o 1 ) = I(o 2 + δ). We thus need to consider the configurations of N 20 whose the central pixel value I(o 2 ) and that of its 4-neighbouring pixel I(o 2 + δ) are different, e.g., I(o 2 ) = 1 and I(o 2 + δ) = 0. Here we set δ = (0, 1), i.e., o 2 + (1, 0) is the right pixel of o 2 . In other words, the pixel values of o 2 and its right pixel o 2 + δ are pre-set. Under such conditions, |C| is reduced to 2 18 . Thanks to the reflection and rotational symmetries, we can again reduce |C| to 124 260 configurations. Then, we use Alg. 1 proposed in Sec. 4.2

Algorithm 2 : 3 B 6 if v ∈ U then 7 UAlgorithm 3 :

 23673 Generation of topologically invariant configuration set P 4 (resp. P 8 ) Input: The DRT graph G p = (V p , E p ) and the set C of 124 260 binary local configurations of size N 20 . Output: The set P 4 (resp. P 8 ). 1 P 4 ← ∅; (resp. P 8 ← ∅;) 2 foreach I C ∈ C do ← T RUE; U ← V p ; S ← {u} where u is the vertex associated to I C in G p ;4 while S ∅ and B = T RUE do 5 Let v ∈ S; S ← S \ {v}; ← U \ {v}; 8 foreach e = (v, w, (p, p ′ )) ∈ E p such that w ∈ U do S ← S ∪ {w}; 9 if ∃ e = (v, w, (o 1 , o 2 )) ∈ E p such that w ∈ U and o 2 is not a simple point in its (4, 8) (resp. (8, 4))-adjacency relations then B ← FALS E; 10 if B = T RUE then P 4 ← P 4 ∪ {I C }; (resp. P 8 ← P 8 ∪ {I C };) Based on Prop. 6 and the sets P 4 and P 8 , the algorithm for characterizing the topologically invariant property of a given binary image I : S → F by a local verification of pixels is given in Alg 3. The algorithm scans I and considers for each pixel p ∈ S its N 20 ( p) with either P 4 or P 8 depending on the binary value of p. Note that I( p) = 1-I( p). Local verification of the topology invariance of a binary image Input: A binary image I : S → F and the sets P 4 and P 8 . Output: Yes if I is topologically invariant and No otherwise. 1 foreach p ∈ S do 2 if I(p) = 1 and I |N 20 ( p) P 4 then return No; 3 if I(p) = 0 and I |N 20 ( p) P 8 then return No; 4 return Yes;

Fig. 7 .

 7 Fig. 7. Some topologically invariant configurations of P 4 (a) and P 8 (b).

Fig. 8 .

 8 (a-f) Some examples of topology-invariant images. (g-i) Three examples of topologyvariant images (left) with their transformed images (right) exhibiting different topologies from their respective original images.

⋆ The research leading to these results has received funding from the French Agence Nationale de la Recherche (Grant Agreement ANR-2010-BLAN-0205 03).


Yukiko Kenmochi

Phuc Ngo

Hugues Talbot

Nicolas Passat

Efficient neighbourhood computing for discrete rigid transformation graph search

Keywords: image registration, discrete rigid transformation, discrete optimization, DRT graph

Rigid transformations are involved in a wide variety of image processing applications, including image registration. In this context, we recently proposed to deal with the associated optimization problem from a purely discrete point of view, using the notion of discrete rigid transformation (DRT) graph. In particular, a local search scheme within the DRT graph to compute a locally optimal solution without any numerical approximation was formerly proposed. In this article, we extend this study, with the purpose to reduce the algorithmic complexity of the proposed optimization scheme. To this end, we propose a novel algorithmic framework for just-in-time computation of sub-graphs of interest within the DRT graph. Experimental results illustrate the potential usefulness of our approach for image registration.

Introduction 1.Discrete rotations and discrete rigid transformations

In continuous spaces (i.e., R n ), rotations are some of the simplest geometric transformations. However, in the discrete spaces (i.e., Z n ), their analogues, namely discrete rotations, are more complex. The induced challenges are not simply due to highdimensionality: indeed, even in Z 2 , discrete rotations raise many difficulties, deriving mainly from their non-necessary bijectivity [START_REF] Nouvel | Characterization of bijective discretized rotations[END_REF]. In this context, discrete rotations -and the closely related discrete rigid tansformations -have been widely investigated.

From a combinatorial point of view, discrete rotations have been carefully studied [START_REF] Nouvel | Configurations induced by discrete rotations: Periodicity and quasiperiodicity properties[END_REF][START_REF] Berthé | Discrete rotations and symbolic dynamics[END_REF][START_REF] Nouvel | Self-similar discrete rotation configurations and interlaced Sturmian words[END_REF], in particular to shed light on remarkable configurations induced by the periodicity of rotations with respect to the discrete grid. At the frontier between combinatorics and algorithmics, the problem of 2D pattern matching under discrete rotations has also been explored [START_REF] Jacob | On discrete rotations[END_REF][START_REF] Amir | Faster two-dimensional pattern matching with rotations[END_REF].

From an algorithmic point of view, efforts have been devoted to effectively compute discrete rotations. In particular, the quasi-shear rotations [START_REF] Reveillès | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF][START_REF] Andres | The quasi-shear rotation[END_REF] were introduced to preserve bijectivity, by decomposing rotations into successive quasi-shears.

Finally, from an applicative point of view, discrete rotations have been used for image/signal processing purposes [START_REF] Richman | Understanding discrete rotations[END_REF][START_REF] Andres | Discrete rotation for directional orthogonal wavelet packets[END_REF]. Other strategies have also been proposed to pre-process 2D images in order to guarantee the preservation of topological properties under discrete rigid transformations [START_REF] Ngo | Topology-preserving rigid transformation of 2D digital images[END_REF].

Recently, we proposed a new paradigm to deal with discrete rotations, and more generally rigid transformations. This paradigm relies on a combinatorial structure, called discrete rigid transformation graph (DRT graph, for short) [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF]. This structure describes the quantification of the parameter space of rigid transformations, in the framework of hinge angles, pioneered in [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF][START_REF] Nouvel | Incremental and transitive discrete rotations[END_REF][START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF].

The DRT graph has already allowed us to contribute to the state of the art on rigid transformations from a combinatorial point of view, by establishing the complexity of "free" [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF] and "constrained" [START_REF] Ngo | On 2D constrained discrete rigid transformations[END_REF] discrete rigid transformations. From an algorithmic point of view, it has been used to characterise topological defects in transformed images [START_REF] Ngo | Topology-preserving conditions for 2D digital images under rigid transformations[END_REF]. Finally, we recently started to explore the applicative possibilities offered by the DRT graph. In particular, we have considered its potential usefulness in the context of image registration [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF].

Registration issues

In the context of image processing, geometric transformations are often considered for registration purposes [START_REF] Zitová | Image registration methods: A survey[END_REF]. Registration is indeed a complex, often ill-posed problem, that consists of defining the transformation that is required to correctly map a source image onto a target image.

Registration is mandatory in various application fields, from remote sensing [START_REF] Schowengerdt | Remote Sensing: Models and Methods for Image Processing[END_REF] to medical imaging [START_REF] Noblet | Recalage d'images médicales[END_REF]. According to the specificities of these fields, registration can implicate different types of images (2D, 3D) and transformations, both rigid and nonrigid. However, the problem remains almost the same in all applications. Given two images A and B, we aim at finding a transformation T * within a given transformation space T. This transformation minimizes a given distance d between the image A and the transformed image T (B) of the image B by T , i.e.

T * = arg min T ∈T d(A, T (B)) (1)
In recent works [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF], we investigated how to use the DRT graph in order to solve this problem in the case of rigid registration of 2D images. The novelty of this approach, with respect to the state of the art, was to provide exact transformation fields, so as to avoid any interpolation process and numerical approximations.

In this context, a preliminary algorithm was proposed for computing a local minimum for Eq. ( 1), thus providing a solution in a neighbourhood of depth k ≥ 1, to the above registration problem. This algorithm strongly relies on the DRT graph, and consists of exploring a sub-graph defined around a given vertex, modeling an initial transformation. Its time complexity was O(m k N 2 ), which is linear with respect to the image size, but exponential with respect to the neighbourhood depth (with m the size of the 1-depth neighbourhood).

Contribution

We propose an improved algorithm (Sec. 3), which dramatically reduces the exponential complexity of that developed in [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF]. Indeed, we show that the k-depth neighbourhood of a DRT graph can be computed with a time complexity O(kN 2 ) (Sec. 4). Experiments emphasise the methodological interest of the proposed approach (Sec. 5).

Introduction to discrete rigid transformation graphs 2.1 Rigid transformation space

In the continuous space R 2 , a rigid transformation is a bijection T : R 2 → R 2 , defined, for any x = (x, y) ∈ R 2 , by

T (x) = cos θ -sin θ sin θ cos θ x y + a 1 a 2 (2)
where a 1 , a 2 ∈ R and θ ∈ [0, 2π[ (T is sometimes noted T a 1 a 2 θ ). In order to apply such rigid transformations on Z 2 , a post-processing digitization is required. More precisely, a digitized rigid transformation T : Z 2 → Z 2 is defined as T = D • T where D : R 2 → Z 2 is a rounding function. In other words, for any p = (p, q) ∈ Z 2 , we have

T ( p) = p q = D • T ( p) = [p cos θ -q sin θ + a 1 ] [p sin θ + q cos θ + a 2 ] (3)
The use of the rounding function D implies that digitized rigid transformations are not continuous within the 3D parameter space induced by a 1 , a 2 and θ. The transformations leading to such discontinuities are called critical transformations. In the space (a 1 , a 2 , θ), the subspace of critical transformations is composed of 2D surfaces Φ pqp and Ψ pqq , analytically defined, for any p = (p, q) ∈ Z 2 and any vertical (resp. horizontal) pixel boundary x = p + 1 2 (resp. y = q + 1 2 ) with p ∈ Z (resp. q ∈ Z), by

Φ pqp : p cos θ -q sin θ + a 1 = p + 1 2 (4)
Ψ pqq : p sin θ + q cos θ + a 2 = q + 1 2 (5)
For a given triplet (p, q, p ) (resp. (p, q, q )), Φ pqp (resp. Ψ pqq ) is called a vertical (resp. horizontal) tipping surface in the parameter space (a 1 , a 2 , θ), and a vertical (resp. horizontal) tipping curve in the 2D plane (a 1 , θ) (resp. (a 2 , θ)).

For an image of size N × N, Φ pqp and Ψ pqq verify p, q ∈ [[0, N -1]] and p , q ∈ [[0, N]]. Examples of tipping surfaces and curves are illustrated in Fig. 1.

Discrete rigid transformation graph

A set of tipping surfaces induces a subdivision of the (a 1 , a 2 , θ) space into classes, each consisting of transformations T a 1 a 2 θ such that (a 1 , a 2 , θ) → T = D • T a 1 a 2 θ is constant. These classes -called discrete rigid transformations (DRTs) -indeed form 3D cells, bounded by tipping surfaces that correspond to discontinuities. By mapping each cell onto a vertex, and each tipping surface piece onto an edge, in a Voronoi/Delaunay paradigm, we can model this subdivided parameter space as a graph, called the DRT graph, as illustrated in Fig. 2.

Definition 1 ([12]

). A DRT graph G = (V, E) is defined such that (i) each vertex v ∈ V models a DRT; and (ii) each labelled edge e = (v, w, f ) ∈ E, where f is either Φ pqp or Ψ pqq , connects two vertices v, w ∈ V sharing the tipping surface f as boundary.

For a given image I, each vertex is associated with a unique transformed image, induced by the DRT corresponding to the vertex. The existence of an edge between two vertices indicates a "neighbouring" relation between the two associated DRTs. More precisely the two transformed images differ by at most one over the N 2 pixels of I; the edge label f provides this information. This allows us to use the DRT graph to produce all the transformed images via successive elementary (i.e., single-pixel) modifications.

Discrete rigid transformation graph and image registration

The registration problem formalised in Eq. ( 1) consists of finding the transformation that best maps a source image onto a target image, with respect to a given distance.

In the discrete framework, the number of transformations is actually finite. In particular, in the case of rigid registration, the solution(s) to Eq. ( 1) can be found within the DRTs exhaustively modeled by the DRT graph. In other words, by considering a brute-force search, a solution, i.e., a global optimum, can be determined. However, the DRT graph G of an image of size N × N, has a high space complexity O(N 9 ) [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF] that induces the same time complexity both for its construction and exhaustive search.

This limits exploration of the whole structure to relatively small images. Nevertheless, as already discussed in [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF], it is possible to perform a local search on G in order to determine a local optimum.

Local search on a discrete rigid transformation graph

To find such an optimum, a local search begins at a given transformation, i.e., a chosen vertex v of G. Then, it moves towards a better solution in its neighbourhood -following a gradient descent -as long as an improved solution can be found. Beyond the choice of the initial vertex -often guided by the application context -the most critical issue is the choice of a "good" search area around this vertex, i.e., a depth of its neighbourhood. In particular, the trade-off is time efficiency versus exhaustiveness.

The neighbourhood of depth 1, noted N 1 (v), actually corresponds to the set N(v) of vertices adjacent to v in G. More generally, neighbourhoods of depth k ≥ 1, also called k-neighbourhoods, are then recursively obtained as

N k (v) = N k-1 (v) ∪ u∈N k-1 (v) N(u) (6)
where N 0 (v) = {v}. Our initial algorithm [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF] was directly mapped on this recursive definition. As a consequence, this approach led to a high time complexity O(m k N 2 ), that is exponential with respect to the depth k of the neighbourhood with vertex degree m, which is supposed to be constant in average (Sec. 4.2). In the next section, we propose a more efficient algorithm, that removes this exponential cost.

k-Neighbourhood construction algorithm

We now propose an algorithm that efficiently computes the part of a DRT graph that models the neighbourhood of depth k around a given vertex. To this end, we need to handle the analytical representation of the cells associated to the DRT graph vertices, inside the subdivided parameter space of (a 1 , a 2 , θ) (Sec. 3.1). Then, we develop a construction strategy that relies on a sweeping plane technique introduced in [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF] (Sec. 3.2). The final algorithm is described and formalized in Sec. 3.3.

Tipping surfaces associated to a discrete rigid transformation

A vertex v of a DRT graph G corresponds to one discrete rigid transformation, that induces a unique transformed image I v obtained by applying this transformation on an initial image I. In other words, for each pixel (p i , q i ) of I v , we know which pixel (p i , q i ) of I transfers its value to (p i , q i ). This correspondence is modeled by the following inequalities deriving from Eq. ( 3)

p i - 1 2 < p i cos θ -q i sin θ + a 1 < p i + 1 2 (7)
q i - 1 2 < p i sin θ + q i cos θ + a 2 < q i + 1 2 (8)
For an image I of size N × N, each of the N 2 pixels generates 4 such inequalities. In the parameter space of (a 1 , a 2 , θ), the obtained 4N 2 inequalities then define a 3D cell, denoted by R v , which gathers all the parameter triplets associated to the discrete rigid transformation corresponding to the vertex v.

When interpreting these inequalities in terms of tipping surfaces/curves (see Eqs. (4-5)), it appears that for each pixel of I v , Eqs. (7-8) 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 )). We derive from Eqs. [START_REF] Reveillès | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF][START_REF] Andres | The quasi-shear rotation[END_REF], that any cell R v is directionally convex along the aaxes [START_REF] Ngo | On 2D constrained discrete rigid transformations[END_REF]. This implies that for any θ value where it is defined, R v is bounded by at least one upper (resp. lower) tipping surface, which constitutes the upper (resp. lower) part of its boundary in each a -direction. This property can be used for constructing a DRT graph locally, or for obtaining topological information from a DRT graph such as a neighbourhood. One may notice that it is sufficient to consider only tipping surfaces of (S + (I v ) ∪ S -(I v )) in order to obtain the k-neighbourhood of v, if k < N.

Sweeping plane algorithm for DRT sub-graph construction

In our new algorithm, the purpose is to build a k-neighbourhood "similarly" to the construction of a 1-neighbourhood in our previous version [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF], that is by using a sweeping plane technique from one value θ v within R v , to both the left-hand and right-hand sides along the θ-axis in the space (a 1 , a 2 , θ).

The differences between this new algorithm and the former are twofold. On the one hand, the range of the considered θ values is wider. Indeed, the sweep must be carried out inside R v but also outside. On the other hand, a larger number of tipping surfaces are considered around R v , while only immediate neighbours were previously involved.

To ease the understanding of this algorithm, we first recall the general idea of the sweeping plane technique.

Given a set S of s 1 vertical and s 2 horizontal tipping surfaces, we aim to construct the DRT sub-graph G corresponding to a given range [θ start , θ end ]. By comparison to [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF], the plane is then swept from θ start to θ end , instead of 0 to 2π. From the very definition of tipping surfaces, this plane is subdivided into (s 1 + 1) × (s 2 + 1) 2D rectangular cells, generated by its intersection with the tipping surfaces of S. More precisely, we have (s + 1) divisions in each a -direction, except at the intersection of tipping surfaces, where a rectangle disappears while a new appears. By observing these rectangle updates during the plane sweeping from θ start to θ end , we can construct the DRT subgraph, where each rectangle corresponds to a vertex while each tipping surface between two rectangles corresponds to an edge. In other words, at each intersection of tipping surfaces, s new vertices and their associated (3s +2) edges are generated, as illustrated in Fig. 3. (The reader is referred to [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF] for more details.)

Our modified algorithm consists of using a topological sweep [START_REF] Edelsbrunner | Topologically sweeping an arrangement[END_REF] in order to find the next closest intersection of tipping surfaces for updating the planar division. We consider at most |S| -2 intersections at each update, by considering only the intersections of consecutive tipping surfaces in their ordered structure in the sweeping plane along each a -axis, and find the closest one among them. After each update, the modifications of such intersections can be performed in constant time. We can also ignore the intersections that are not in the range between θ start and θ end . In particular, since we have |θ end -θ start | 2π, the number of intersections can be considered a small constant.

Hereafter, we denote this specific procedure by S weep(S, θ start , θ end ), and we write G = S weep(S, θ start , θ end ) for the output DRT sub-graph.

k-Neighbourhood construction

Finding the neighbouring vertices and edges of a given vertex v with depth k, is actually equivalent to constructing the DRT sub-graph containing those vertices and edges around v. Here, we assume to know a value θ v lying into R v , and we use it as initial value of the sweeping algorithm. The plane is thus swept twice, in the space (a 1 , a 2 , θ), along the two directions of the θ-axis.

The key-point is how to limit the construction of the DRT sub-graph. For this purpose we verify, for each generated vertex u, its neighbourhood depth t v (u) with respect to v. If t v (u) > k for all vertices in the current sweeping plane, the process ends.

When a vertex u is created, its depth t v (u) depends on that of the two vertices w 1 and w 2 to which it is adjacent in the a -direction of the tipping surface intersection. We then have t v (u) = 1 + min{t v (w 1 ), t v (w 2 )}. (An iterative backtracking process is also needed to update the depth of w and its successive neighbours, whenever t v (w ) > t v (u) + 1.) 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 .

4

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

5 Initialize V = ∅, E = ∅ 6 Initialize θ prev = θ v 7 repeat 8 for = 1, 2 do 9
Find the next left intersection θ + of f + (resp. θ -of f -) with the other surface in S for θ < θ prev .

10

θ next = min{θ + 1 , θ - 1 , θ + 2 , θ - 2 } -with 1 11 ∆F = S weep(S 1 ∪ S 2 , θ prev , θ next ) 12 if ∃u ∈ ∆V, t v (u) ≤ k then 13 F = F ∪ ∆F, θ prev = θ next 14 if the next intersecting surface with f + (or f -) is in S then 15
Exchange their order in S . Replace f + (or f -) in S with the new intersecting surface.

18 until ∀u ∈ ∆V, t v (u) > k;
In each a -direction, by considering the (k + 1) closest tipping surfaces around R v , we can obtain all the vertices u such that t v (u) ≤ k. In the θ-direction, we need to check if t v (u) > k for all vertices u in the current sweeping plane; if so, the sweeping ends. The global process is described in Alg. 1. (Note that the algorithm describes only the k-neighbourhood construction in the left-hand side along the θ-axis, but the righthand side can be constructed similarly.) The first loop (Lines 1-4) initializes the set of tipping surfaces that are needed to generate the k-neighbours of a given DRT v. We obtain 2(k + 1) vertical (resp. horizontal) tipping surfaces close to R v at θ = θ v , and sort and store them in the lists S . In the second loop (Line 7), we first verify how long we can keep the same tipping surface sets S (Lines 9-10), and then build a DRT sub-graph by using the S weep algorithm for this verified θ interval (Line 11). After verifying if there still exists a generated vertex whose neighbourhood depth is ≤ k (Line 12), we update the tipping surface sets S for the next interval (Lines 14-17).

Obviously, F is not the smallest sub-graph G including the k-neighbours of v. To obtain G from F, we simply keep vertices whose neighbourhood depth is ≤ k. In order to obtain the initial 2(k + 1) vertical (resp. horizontal) tipping surfaces of S , the time complexity is O(N 2 ) for Line 2; O(N 2 ) for Line 3 on the average case 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 4, respectively. Then, we carry out the plane sweep for each updated S 1 ∪ S 2 . For each iteration in the loop, the most costly parts are Lines 9 and 11, which require O(N 2 ) and O(k 2 ), respectively. The next question concerns the number of updates for S 1 ∪ S 2 . If m is the degree of any vertex u of a DRT graph, this update number can be estimated as m(2k + 1), since the union of R u for all u in the k-neighbourhood of a given vertex v contains at most 2k + 1 adjacent R u in the θ-direction. Therefore, the time complexity is O(mkN 2 ) for this iterative plane sweep loop.

The time complexity of Alg. 1 is thus O(mkN 2 ), which is significantly lower than that of our previous algorithm [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF], namely O(m k N 2 ). We observe, in the next section, that m can be estimated as a low constant value, leading to a final complexity of O(kN 2 ).

Average degree of DRT graphs

The DRT graph space complexity for an image of size N × N is O(N 9 ), both for vertices and edges [START_REF] Ngo | Combinatorial structure of rigid transformations in 2D digital images[END_REF]. In other words, the number of vertices and that of edges grow at the same rate. We can then infer that m is actually bounded, independently of N.

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. Indeed, 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 v-e+ f = 2, where v, e and f are the number of (0D) vertices, (1D) edges and induced (2D) cells, respectively. From the very definition of the DRT graph, 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. This is confirmed by the experimental analysis, illustrated in Fig. 4(a).

In practice, the maximal degree of the vertices within a DRT graph also remains close to this average value. Indeed, the histograms depicted in Fig. 4(b) show that the 2D DRT vertex degrees converge rapidly to a stable distribution that contains mainly degrees of value 3 and 4 (with a maximal value experimentally identified at 8). More qualitatively, Fig. 5 illustrates the distribution of these degrees of a 2D DRT graph.

Experiments

Iterated local search was applied to image registration. In this section we validate Alg. 1 in practice, and we observe its local behaviour when varying k. For simplicity, we use the same experimental settings as in [START_REF] Ngo | Discrete rigid transformation graph search for 2D image registration[END_REF], i.e.two input binary images and a signed distance [START_REF] Boykov | An integral solution to surface evolution PDEs via geo-cuts[END_REF] for Eq. [START_REF] Nouvel | Characterization of bijective discretized rotations[END_REF]. In order to find an initial transformation, we use the first and second order central moments of a binary shape [START_REF] Flusser | Moments and Moment Invariants in Pattern Recognition[END_REF]. Experiments are carried out with different neighbourhood sizes, k = 1, 3, 5, 10, 15 on binary images of size 53 × 53 from the initial transformation, as illustrated in Fig. 6. We can observe in Fig. 7 that the locally optimal distance improves when we use a larger neighborhood, which is coherent in a gradient descent paradigm.

Conclusion

We have significantly improved the time complexity of the process of computing a neighbourhood of given depth within a DRT graph, without requiring the computation of the whole graph. This time complexity may be reduced in some cases, in particular if the image is binary by dealing only with the pixels that compose the binary object border. The proposed applications only validate our approach as a proof of concept. Nevertheless, an exacti.e., numerical error-free -strategy is novel in the field of image registration and may open the way to new image processing paradigms. In future work we will explore the notion of DRT graph in Z 3 .

Fig. 1 .

 1 Fig. 1. (a) Tipping surfaces in the space (a 1 , a 2 , θ), and (b) their tipping curves [16].

Fig. 2 .

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

Fig. 3 .

 3 Fig. 3. 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 DRT graph in each plane (in green: new vertices and edges in the second and third planes).

Algorithm 1 :

 1 k-neighbourhood computation (in the left-hand side along the θ-axis) Input: A DRT v (or its associated image I v ); a positive integer k. Output: The DRT sub-graph F = (V, E) containing the k-neighbours of v. 1 for = 1, 2 do 2 Determine the tipping surfaces associated to v: S + (I v ), S -(I v ) (Sec. 3.1).

3

 3

Fig. 4 .

 4 Fig. 4. (a) Average vertex degree in a 2D DRT graph. (b) Normalised vertex degree distribution in a 2D DRT graph.

Fig. 5 .

 5 Fig.5. 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).

Fig. 6 .

 6 Fig. 6. Input images and results of the iterated local search for image registration. (a) Reference image, (b) target image, and (c) the initial transformed image of (b) 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 respectively. Note that in (c-h), pixels are coloured if they are different from those in (a); yellow (resp. red) pixels are white (resp. black) in (c-h) and black (resp. white) in (a).

Fig. 7 .

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

The research leading to these results has received funding from the French Agence Nationale de la Recherche (Grant Agreement ANR-2010-BLAN-0205).