J S Bailly
email: bailly@teledetection.fr

F Levavasseur
email: levavasseur@supagro.inra.fr

P Lagacherie
email: lagacherie@supagro.inra.fr

A spatial stochastic algorithm to reconstruct artificial drainage networks from incomplete network delineations

Keywords: channels, ditches, mapping, graphs, simulated annealing, simulation, uncertainties, remote sensing, random walks

A spatial stochastic algorithm that aims to reconstruct an entire artificial drainage network of a cultivated landscape from disconnected reaches of the network is proposed here. This algorithm uses random network initialisation and a simulated annealing algorithm, both of which are based on random pruning or branching processes, to converge the multi-objective properties of the networks; the reconstructed networks are directed tree graphs, conform to a given cumulative length and maximise the proportion of reconnected reaches. This algorithm runs within a directed plot boundaries lattice, with the direction governed by elevation. The proposed algorithm was applied to a 2.6-km 2 catchment of a Languedocian vineyard in the south of France.

The 24-km-long reconstructed networks maximised the reconnection of the reaches obtained either from a hydrographic database or remote sensing data processing. The distribution of the reconstructed networks compared to the actual networks was determined using specific topographical and topological metrics on the networks. The results show that adding data on disconnected

Introduction

Artificial drainage networks of agrarian landscapes include connected linear features, such as channels, culverts, roads and reshaped gullies, that result from the constant efforts of farmers to adapt landscapes to the constraints of agricultural production. They constrain water flow paths in landscapes, which greatly impacts the environmental functioning of agrarian areas: drainage networks can be ecological hot-spots and corridors [START_REF] Carpentier | Trends of a bitterling (rhodeus sericeus) population in a manmade ditch network[END_REF][START_REF] Croxton | Evaluation of automatic road extraction[END_REF][START_REF] Pita | Conserving the cabrera vole, microtus cabrerae, in intensively used mediterranean landscapes[END_REF] and control many fluxes, especially of water and pollutants [START_REF] Bouldin | Vegetative and structural characteristics of agricultural drainages in the mississippi delta landscapes[END_REF]Dages et al., 2008).

When scaling up from an agricultural plot scale, the structure of drainage networks appears to be highly variable in space [START_REF] Lagacherie | An indicator approach for describing the spatial variability of human-made stream network in regard with herbicide pollution in cultivated watersheds[END_REF].

Therefore, most environmental diagnoses and impact assessments at the catchment scale in agrarian areas now take these networks and their spatial variabilities into account through physical modelling [START_REF] Lagacherie | Geomhydas: A landscape discretization tool for distributed hydrological modeling of cultivated areas[END_REF][START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF][START_REF] Varado | Contribution au développement d'une modélisation hydrologique distribuée. application au bassin versant de la donga au bénin[END_REF]. However, to date, these approaches have only been used for small areas where major investments were made in terrain surveys to map drainage networks [START_REF] Lagacherie | An indicator approach for describing the spatial variability of human-made stream network in regard with herbicide pollution in cultivated watersheds[END_REF], which limits the generalisation of these studies. Consequently, the extension of environmental agrarian landscape diagnoses to wider areas with a reasonable cost is largely dependent on our ability to map such linear elements.

Similar to most surface linear mapping methods (Amini et al., 2002), remote sensing represents a candidate method for artificial drainage network mapping [START_REF] Vassilopoulos | Geoinformation Technologies for Geo-Cultural Landscapes: European Perspectives[END_REF]. Methods allowing the delineation of drainage networks in natural areas from regular grid digital terrain models (DTMs) have existed since 1984 [START_REF] O'callaghan | The extraction of drainage networks from digital elevation data[END_REF]; these methods use drainage algorithms to search for the flow direction in each pixel [START_REF] Wilson | Comparison of the performance of flow-routing algorithms used in gis-based hydrologic analysis[END_REF] and build aggregated drainage areas accordingly (e.g., vectorised thresholded drainage areas usually form the final drainage networks) [START_REF] Soille | An efficient algorithm for drainage network extraction on dems[END_REF].

For cultivated landscapes, drainage algorithms used on DTMs are unable to represent anthropogenically modified overland flow paths (Duke et al., 2006;Garcia and Camarasa, 1999), even when using high resolution airborne scanner laser DTMs [START_REF] Comoretto | Analyse des chemins de l'eau d'un petit bassin rural sur mnt thr[END_REF].

Instead, methods delineating network reaches from images extracted from whole networks have previously been assessed [START_REF] Bailly | Agrarian landscapes linear features detection from lidar: application to artificial drainage networks[END_REF][START_REF] Vassilopoulos | Geoinformation Technologies for Geo-Cultural Landscapes: European Perspectives[END_REF]. However, as for most raw image processing results, confusion in remote sensing detection provides an incomplete view of a network by producing sets of unconnected and uncertain reaches.

Consequently, methods that aim to reconstruct an entire likely connected drainage network from an incomplete delineated network must be proposed.

These methods could be used to reconnect reaches detected by remote sensing, as well as to extend known downstream regions of networks existing in geographical databases. In the case of mapping explicit uncertainties, these methods might allow assessment of the ultimate utility of specific remote sensing processes for network mapping or the required amount of known downstream reaches of networks in geographical databases. Furthermore, with regard to the use of drainage network maps for modelling the geometry of fluxes [START_REF] Beven | The future of distributed models: Model calibration and uncertainty prediction[END_REF][START_REF] Beven | Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the glue methodology[END_REF], methods are also needed that represent uncertainties in network reconstruction, i.e., uncertainties in network mapping.

Methods that aim to reconstruct line networks (e.g., roads and hydrographic networks) from incomplete network delineations from remote sensing [START_REF] Stoica | A gibbs point process for road extraction in remotely sensed images[END_REF][START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF] or vectorised cartographic elements [START_REF] Baltsavias | Automated updating of road databases from aerial images[END_REF][START_REF] Mariani | Image Analysis Applications and Computer Graphics[END_REF] have previously been proposed. The former method uses a random polyline stochastic simulation process derived from the Candy model (Descombes et al., 2001), and the latter uses a deterministic minimisation process [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF][START_REF] Prim | Shortest connection networks and some generalizations[END_REF]. Both reconnect networks with geometric rules, favouring long, slighty curved polylines without overlapping or aligned segments. However, at the end of the process, both methods provide a single estimated network, without computation of uncertainties.

Considering this last point, we chose first to develop a reconstruction algorithm to: 1) randomly reconstruct various possible artificial drainage networks from an incomplete delineation of the network and 2) reconstruct networks taking into account a number of global network geometrical and topological properties. The algorithm employed here is a spatially simulated annealing algorithm that was chosen for its ability both to better explore all network possibilities and to escape from local minima during convergence [START_REF] Vangroenigen | Constrained optimization of spatial sampling using continuous simulated annealing[END_REF]. This algorithm generates networks corresponding to directed trees (Fairfield and Leymarie, 1991) and reconnecting reaches of the incomplete delineated network. The reconstructed networks correspond to connected sub-graphs of the directed lattices of agricultural plot boundaries. The proposed algorithm falls within the scope of spatial conditional simulation methods [START_REF] Lantuejoul | Geostatistical Simulation: Models and Algorithms[END_REF]. The sampled reaches of the network conditioning the simulation can come from ground observations, remote sensing detection or existing geographical databases. The direction of the graph comes from elevation gradients. Uncertainties in mapping networks are simulated through a set of equi-probable generated networks.

To stringently evaluate how these uncertainties impact network morphology in a catchment space, methods measuring how similar the reconstructed networks are to the actual networks need to be employed. In network mapping, this is usually carried out using only geographical overlap criteria [START_REF] Molly | Automatic mapping of valley networks on mars[END_REF][START_REF] Croxton | Evaluation of automatic road extraction[END_REF]. However, this method does not permit the assessment of how realistic the morphology of the reconstructed network is. Thus, specific complementary similarity metrics of networks were developed for this purpose. This paper is organised into four sections: 1) a description of the developed methods, i.e., the reconstruction algorithm (its assumptions, principles and general parameters) and the network topological and and topographical dissimilarity metrics are first exposed; 2) a case study on the drainage ditch network of a 2.6 -km 2 vineyard catchment located in the south of France, with different sets of data conditioning the network reconstruction and two simulated annealing convergence criteria being used for network reconstruction is secondly presented; 3) comparisons of topological and and topographical dissimilarities between the actual network and some reconstructed networks are performed; and 4) we conclude and expand on the limits of the proposed method to properly assess the implications of the findings produced from an incomplete delineation network from remote sensing. denotes n i as the start-vertice and n - i denotes n i as the end-vertice.

Materials and methods

The reconstruction algorithm

Assumptions

We first considered that the graph d that we want to simulate is a connected sub-graph of the valued, connected, planar lattice representing the cultivated plot boundaries within a cultivated landscape. ADN reaches are, in fact, almost always located on plot boundaries. This lattice is denoted similarly as l = (x l , u l ).

l is a directed graph (digraph), i.e., a directed lattice. The direction is given by the elevation gradient between edge vertices; it corresponds to the upstream-downstream direction. However, considering that elevation is uncertain due to the noise σ z parameter, an ambiguous direction for a given plot boundary ending on the two vertices n i , n j can be obtained. In this case, we chose to retain the two edge directions, which meant that two distinct edges, v and v , were considered possible between vertices n i , n j :

∀n i , n j ∈ x 2 l / |zn i -zn j |<σz , v = (n + i , n - j ) ∈ u l , v = (n + j , n - i ) ∈ u l (1)
As has been done in most ADN hydrology studies [START_REF] Lagacherie | Geomhydas: A landscape discretization tool for distributed hydrological modeling of cultivated areas[END_REF][START_REF] Moussa | Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events[END_REF][START_REF] Varado | Contribution au développement d'une modélisation hydrologique distribuée. application au bassin versant de la donga au bénin[END_REF], we considered that the d graph is also a directed tree (ditree) with a root r corresponding to the physical outlet.

To assimilate any type of existing data for ADN locations, we also consider that an unconnected sub-graph exists:

f = (x f , u f ) ∈ d. u f denotes
a set of fixed edges. This can be produced from an existing geographical database related to a hydrographic network for downstream reaches, from a remote sensing detection process, or from any field survey information. As an extremal case, f can only contain the r vertex, the fixed outlet. Additionally, the subsets of fixed edges u f 1 and u f 2 , such as

u f 1 + u f 2 = u f , are
also distinguished. u f 2 denotes the set of fixed edges with a fixed direction (typically obtained from a geographical database), regardless of the slope gradient. u f 1 denotes the set of fixed edges with an uncertain direction (typically obtained from a remote sensing detection process).

Spatial simulated annealing algorithm general scheme

Simulated annealing is carried out for the optimisation of stochastic algorithms based on the Metropolis-Hastings (MH) heuristic. Basically, it allows the optimisation of paths through costly points when the temperature is still high in the first iterations of the algorithm. It is an heuristic that is usually employed for spatial problems, and it permits a better exploration of space from an initial solution [START_REF] Vangroenigen | Constrained optimization of spatial sampling using continuous simulated annealing[END_REF].

The proposed formulation of the algorithm can be written in pseudo-code as follows:

In this algorithm, the parameters are the generally employed simulated annealing parameters: k max is the maximum number of iterations; e max is Algorithm 1 Reconstruction digraph algorithm 1: procedure Simulated.Annealing(d 0 , e max , k max , T ini , g)

2: k ← 0 3: d ← d 0 initial state 4:
e ← energy(d 0 ) initial energy 5: In the following, the initialisation d 0 , the production of a random neighbour state neighb, the decreasing temperature law and the definition of the energy function are successively detailed. Prior to this, the random walk algorithm, on which the two former algorithms are based, is briefly presented.

T ← T ini initial temperature 6: while k < kmax & e >
else if exp ( e-enew T ) > random(U [0, 1]) then MH heuristic 12: d ← d new ; e ← e new

Random walk

Here, a random walk process consists of random sequential and independent paths in a directed lattice starting from an upstream vertex and moving down to a set of downstream end vertices [START_REF] Spitzer | Principles of the Random Walk[END_REF].

Random branching initialisation

A random initial state for a digraph d 0 is obtained using a sequential branching process based on random walks within lattice l. The process begins from the single graph r. At each iteration based on ordered x f vertices related to ascending elevation, the process connects the x f considered vertex to the digraph obtained in the previous iteration by a random walk. The digraph grows to obtain a given ratio of u f edges connected to the root r. In this process, an edge belonging to the (u l -u f 2 ) set is first randomly selected (selection among the uncertain reaches), and two cases are therefore distinguished:

• If this selected edge already belongs to the u d digraph set of edges, then the directed sub-graph (branch) upstream of the ending vertex of this edge is first pruned from digraph d to form d new . Second, a re-branching process is selected through a binary process with a probability of 0.5. If selected, this re-branching process connects each of the u f edges belonging to the upstream sub-graph with same sequential branching process as was used for the initialisation of d 0 , and this new branch is added to d n ew.

• If the selected edge does not belong to (u d ), a random walk starting from its end node, reaching one of the d vertices, is performed.

Consequently, a new branch is added to form d new .

Due to its construction, this neighbour process allows the exploration of the entire l lattice for the reconstruction of d.

Energy functions

The choice of the energy function in the algorithm1 governs the objective properties for the reconstructed networks. Here, we chose a bi-objective energy function. Our aim was that a reconstructed network should:

• converge to a target overall cumulative network length related to the extent of the catchment: this defines the catchment objective catch ; by using this convergence criterion, it is assumed that in future cases, the drainage network density or the cumulative length of a catchment can be predicted from environmental or anthropogenic variables (e.g., soil, catchment morphology and cultivated plot areas. . . ).

• reconnect a maximum number of fixed edges at the extent of the catchment: this defines the reconnection ratio objective ratio, equal to the fraction between the number of reconnected fixed edges and the number of fixed edges;

The multi-objective energy function is energy = max(e 1 , e 2 ), where e 1 , e 2 denote component energy objective functions:

catchment energy: e 1 = |L d -catch| σ 1 with σ 1 = catch * γ 1 reconnection ratio energy: e 2 = |L f ∈d -L f | σ 2 with σ 2 = L f * γ 2
In these formulations, Due to the construction of these formulations, all of the energies have normalised magnitude ranging from 0 to 1.

L d = v∈u d length(v), L f = v∈u f length(v), L f ∈d v∈(u d &u f ) length(v)

Simulated annealing decreasing temperature law

The simulated annealing initial decreasing temperature law follows the commonly employed equation:

T k = T ini * (α k ) (2)
where:

• δ = max(1, (L l -catch) catch * 0.5) • T ini = δ 0.51
The initial temperature parameter δ was fixed to a 'great' cost, divided To assess the accuracy of the proposed ADN mapping method with respect to the f data or algorithm parameters, we computed the ditree network similarities based on scalar metrics. Theses metrics were used to compare each of the simulated ADNs to the actual one and, consequently, the forest dissimilarity distribution was computed.

As the objective function in algorithm 1 already uses global network geometry criteria (the drainage density), we chose to develop two complementary network similarity metrics based on: 1) network topography and 2) network topology.

ADN topographical similarity metric

The topography of an ADN is fully contained in the slope of its network reaches. To summarize these slope distributions, any single synthetic slope statistic (e.g., mean and quantile) performs poorly. Therefore, we chose to use a metric based on the cumulative distribution of edge slopes. For a given network, we calculated the curve representing the cumulative distribution of slopes as follows: 1-for a given slope value θ, the proportion of the network length having slope smaller or equal than θ was computed; 2-this was repeated for discretized and crescent values of the slope ranging from -90 to 90 degrees (see Figure 2).

For two networks, the topographical similarity corresponds to the area between the two curves. This metric h t can be defined as:

h t (d 1 , d ref ) = | 90 -90 F i (θ)dθ - 90 -90 F ref (θ)dθ|, ( 3
)
where θ denotes the slope in degrees and F i (θ) denotes the empirical cumulative distribution function of the network slopes in relative cumulative lengths.

ADN topological similarity metric

Methods for the computation of the similarity between tree topologies have previously been developed in computational science to answer prob-lems related to character chains, phylogenetic tree comparisons or the description of deciduous tree architecture (Ferraro and Godin, 2000). These methods are all based on the Robinson-Foulds metric (Day, 1985) or Levenshtein or Jaro-Winkler distances [START_REF] Jaro | Advances in record linking methodology as applied to the 1985 census of tampa florida[END_REF]. These metrics are dedicated to valued trees with qualitative node labels. However, this is clearly not the problem in our case. In hydrology, topological characterisation of drainage networks is usually based on the Horton-Strahler, Shreve or Tokunaga upstream-downstream edge or node taxonomies [START_REF] Barker | Quantitative morphometry of the branching structure of trees[END_REF].

However, for a given ditree, these taxonomies cannot be easily reduced to a single value: it is usually represented in vectors or matrices as a branching matrix based on bi-order distributions, i.e. distributions of the twostrahler orders for stream joining downstream (Janey, 1992, chapitre3).

To address this problem, we chose to develop the simple similarity metric h c , as proposed initially by Moslonka-Lefebvre et al. ( 2010) for comparison of gully networks in badlands extracted from remote sensing data [START_REF] Thommeret | Extraction of thalweg networks from dtms: application to badlands[END_REF]. This metric is the distance between the Strahler barycenters of networks, where the Strahler barycenter of a network is the weighted barycenter of the network edge upstream nodes, with weights equal to Strahler orders:

h c (d 1 , d 2 ) = distance(G 1 , G 2 ) (4) • x G 1 = 1 i w 1 i i w 1 i x 1 i ; y G 1 = 1 i w 1 i i w 1 i y 1 i are the coordinates of the d 1 Strahler barycenter tree; • x G 2 = 1 j w 2 j j w 2 j x 2 j ; y G 2 = 1 j w 2 j
j w 2 j y 2 j are the coordinates of the d 2 Strahler barycenter tree;

• w 1 i , w 2 j denotes the Strahler order of a node equal to that of the edge for which the node is the upstream for trees d 1 , d 2 respectively. x 1 i , x 2 j and y 1 i , y 2 j denote the coordinates of the edge upstream vertices of trees d 1 , d 2 respectively. This unique metric combines the topology, geometry and geography of the network. It is further simply considered as the topological similarity metric.

Case study

The Muscat catchment and drainage network

Study site

The study area used here is the 2.6-km 2 Muscat catchment located in the Languedoc region (south of France). This area was selected as being representative of the Languedocian vineyard landscape, where linear elements play a major role in hydrological transfers. The area is successively composed of plateaus, sloping areas with terraces, gentle foot slopes formed from alluvial and colluvial deposits and a flat depression from upstream to downstream. The study area includes the Roujan elementary catchment, which has been equipped for conducting measurements of pollution, erosion and flooding since 1992 [START_REF] Voltz | Omere: Observatoire mediterraneen de l'environnement rural et de l'eau: Impact des actions anthropiques sur les transferts de masse dans les hydrosystemes mediterraneens ruraux. proposition d'observatoire de recherche en environnement[END_REF]. The drainage network drains cultivated plots and prevents erosion in the case of major runoff events. The positioning of the reaches that comprise most of this network within the landscape is highly variable and is highly dependent on reach functions. When the reach function is to capture surface waters, as in the case of high slopes, they are perpendicular to the terrain slope, especially at the footstep of terrassettes. When the reach function is water transport, the reaches follow the slope direction. When reach function is to drain subsurface waters, there is no evident relationship with the local relief.

The actual Muscat drainage network

The actual drainage network of the studied catchment was exhaustively surveyed in 2004. It contains 239 edges and 240 nodes with 123 sources and 116 confluences (Fig. 3-a). Its cumulative length is 23,995 m that defines the catch parameter value.

The plot lattice and DTM data

For this catchment, a cultivated connected plot lattice containing 2,645 edges and 1,079 nodes was used (Fig. 3-b). The elevation was assigned to each vertex of this lattice based on a 5-m resolution DTM having a random error of σ z = 1m for the standard deviation. There is no internal sink within this lattice that can trap random directed walks.

The ADN conditioning data

We used the whole set of known edges, f , to condition the reconstruction, as depicted in Figure 3-c. This set was produced from:

• remote sensing detection based on LIDAR raw data (red lines). Remote sensing provides disconnected reaches, which may represent false detection, with a probability of 0.17 of false detection for upstream reaches that were not previously mapped in a geographical database.

The cumulative length of the reaches is 15,400 m, i.e., 47% of the actual network cumulative length. This set of fixed reaches corresponds to the f 1 edge set because the direction of the reaches is unknown.

• the c BD-Topo-geographical database for hydrographic networks for the main downstream reaches (blue). Their cumulative length is 2,500 m, i.e., 7.5% of the actual network. This set of fixed reaches corresponds to the f 2 edge set, which has a known direction.

• r which is the known outlet of the network that we want to reconstruct.

Although it is not appropriate (no edge), we further consider that r can be the extreme case for the f conditioning data: only the outlet is fixed.

Results: reconstructed ADN on the Muscat catchment

Convergence efficiency

Beginning from the cultivated plot lattice shapefiles, algorithms were developed in the R statistical computing program [START_REF] Ihaka | R: A language for data analysis and graphics[END_REF] using the shapefiles and e1076 libraries [START_REF] Team | R: A Language and Environment for Statistical Computing[END_REF]. It took 9, 7.2 and 48 minutes on average to reconstruct a network using the (f = r),

(f = f 2 ) (f = f 1 + f 2 )
sets for conditioning, respectively, on a dual core PC with 1 GB RAM. Additionally, r, f 2 , f 1 + f 2 denote the three exclusive cases of f fixed edges used to condition the reconstruction.

The e max , γ 1 , γ 2 , k max parameters were adapted depending on the f conditioning data:

• when using r or f 2 to condition the reconstruction, the convergence parameters were γ 1 = 1, γ 2 = 1, e max = 0.01 and k max = 500.

• when using f 1 + f 2 to condition the reconstruction, it took a long time to achieve convergence due to both considerable constraints related to lattice space exploration and a multi-objective (e 1 , e 2 ) dilemma [START_REF] Kung | On finding the maxima of a set of vectors[END_REF]. A Pareto optimum providing a compromise between the duration of and the objective of accuracy for the simulation was found with the following convergence parameters: γ 1 = 1, γ 2 = 4 for e max = 0.05 and k max = 2, 000.

In Figure 4, the evolution of energy within the simulated annealing process using f = (f 1 + f 2 ), as well as according to the reconstructed digraph, is depicted for a simulation example at three intermediary states. In Figure 4, the cumulative length of the initial digraph d passing through at least 95 % of the f edges covers most of the lattice l and, consequently, has a high energy e due to its overly long cumulative network length. At the intermediary state (k = 30), the energy is still high (around 0.2) due to a weak f edge reconnection, and the graph d is segmented in spatial clusters due to a great amount of pruning occurring when the temperature remains high enough. For larger iterations (k > 200), the energy slowly decreases with local pruning and branching to converge to the desired maximum level of energy (e max = 0.05) at iteration k = 308.

Reconstructed network variability

Three forests of 100 reconstructed networks were built using r, f 2 , f 1 + f 2 to condition the reconstruction. Some of the networks sampled from these three forests are depicted in Figure 5. In this Figure, the reconstructed networks look realistic compared to the actual network (Figure 3-b). However, differences appear between cases r, f 2 and f 1 + f 2 . The topology seems less realistic for r than for f 1 + f 2 .

However, the reconstructed networks look more realistic in the parts of the catchment where the actual main branches belong to the f 2 edges.

Conditioning data effects on network similarities

Similarities of the reconstructed network forests using the different sets of fixed edges r, f 2 and f 1 + f 2 were computed. The h t , h c similarity distributions can be deduced from the graphs in Figure 6. The top part of Figure 6 shows that the h t topographic similarity distribution is quite similar for r and f 2 . For f 1 + f 2 , the slope distribution appears to be of higher magnitude where the slope is greater than 0.2 rad.

In all cases, the reconstructed network slopes are not totally centred on the actual slope distribution and overestimate very low and negative slopes. The bottom part of Figure 6 shows that the h c topological similarity distribution is once again quite similar for r and f 2 , but with a lower magnitude for f 2 . For f 1 + f 2 , the topological dissimilarities appear to be better centred, despite an apparently higher magnitude.

Mean and standard deviation statistics were computed from these different similarity distributions (Table 4.2) and then tested for statistically significant differences, using a two sample comparison Welch T test (Table 4.2). These statistics confirm that the use of f 2 slightly improves the topological and topographical similarities of the networks. Despite having a higher ratio of fixed reaches compared to f 2 , the use of f 1 + f 2 significantly improves the topological similarity of the networks but a significantly higher topographic dissimilarity also arises (Table 4.2).

f : r f 2 f 1 + f 2 Actual edges
Probability value H 0 :

r vs f 2 r vs f 1 f 2 f 2 vs f 1 f 2
Topographical metric h t 3.05e-02 8.96e-44 5.91e-38

Topological metric h c 2.36e-03 3.03e-15 6.33e-10 This is probably due to the locations of the f 1 edges corresponding to false remote sensing detection [START_REF] Bailly | Agrarian landscapes linear features detection from lidar: application to artificial drainage networks[END_REF], with more being located in the terraced hillslides of the Muscat catchment and favouring higher slopes in the networks. This could also be produced from the necessary e max difference (relaxation), favouring added edges upstream on terraced hillslides.

From the reconstructed forests and similarity distributions, a median network corresponding the reconstructed network with minimized median orders on similarity metrics was computed. This median network is showed using a) r, b) f 2 , c) f 1 + f 2 on the left of figure 7. In addition, uncertainties in network reconstruction was computed through the frequency in reconstructed networks for each one of the cultivated plot boundaries lattice edge. These uncertainties are also mapped using a) r, b) f 2 , c) f 1 + f 2 on the right on figure 7. On this figure, the median tree using f 1 + f 2 looks the more realistic and close to the actual one. As expected, uncertainties are higher on upstream areas corresponding to terraced hillslides.

Conclusions

In this study, a simulated spatial annealing algorithm for the stochastic reconstruction of artificial drainage networks was proposed. This algorithm uses the plot boundary lattice of a cultivated landscape directed by DTM upstream-downstream ditree structure were produced. They follow a given cumulative length and maximise the reconnection of unconnected known reaches. Randomness in the algorithm comes from upstream-downstream random walks along the cultivated plot lattice and according to random pruning or branching processes in the network.

This algorithm can be used to complete or interpolate from known regions of networks produced from ground observations, geographical databases or remote sensing detection. It also permits the representation of uncertainties in drainage network mapping through multiple realisations of the algorithm.

The results of applying the algorithm to the Muscat catchment with or without using edges detected from airborne scanner laser, from the c BD-Topo database, showed quite realistic reconstructed networks in all cases; the simulations are highly constrained by the 3D structure of the plot lattice. However, a quantitative study of the geometric, topographical and topological properties of the reconstructed networks using different criteria than were used as the objectives for simulated annealing resulted in equivocal conclusions. Adding uncertain information from remote sensing of network locations does not always benefit the model, as it increases the biases of some parameters of the network and reduces others, while adding considerable constraints related to convergence and a long duration to the reconstruction process. It appears that the location, rather than the ratio, of the unconnected true reaches is more important. This assumption could be assessed in future studies using the proposed reconstruction algorithm.

Another question regarding the realistic character of the simulations obtained using this algorithm is related to the final usage of the mapped drainage networks. In this case, only descriptive criteria were used, but cri-teria based on the hydrological functioning of the simulated networks also need to be assessed. In other words, do uncertainties in drainage network maps resulting from the available data and the proposed algorithm have a significant influence when they are propagated into a given hydrological model ?

The proposed algorithm is generic and can be extended to other networks, to other catchments and even using objectives related to other network properties. Its main requirements are a cultivated plot lattice that can be valued with elevations and known outlets. It also requires a value of the cumulative length of the network to be simulated. This drainage density is expected to be predicted spatially in future studies.

Acknowledgement(s)

The authors are grateful to the INSU-ECCO-PNRH French national program in hydrology for its general support of this study, which is included in the MOBHYDIC project. They are also grateful to the Languedoc-Roussillon region for the funding of this study through a Ph.D program.

Aarts, E., Korst, J., vanLaarhoven, P.

  An artificial drainage network (ADN) is considered to be a physical network corresponding to a directed and valued graph d with an edge value equal to the reach length of an according physical network and a vertex value equal to the elevation. Using the usual graph theory notations, d = (x d , u d ) where x d denotes the set of d vertices and u d the set of d edges. An edge corresponds to a directed or ordered pair of vertices (n + i , n - j ), where n + i

  the energy objective criteria to reach; and T ini is the initial temperature. g(T, k) denotes a given decreasing law of temperature. d 0 corresponds to the initial state for digraph reconstruction. The random(U [0, 1]) function gives a sample from the uniform distribution. At the end of the algorithm 1, a digraph d connected to root r is obtained satisfying the e max criterion. Due to its construction, few cycles can occur in this digraph. To ensure a final ditree structure for d, the cycles are broken in a final deterministic uncycling process that disconnects the longest paths from the root (in edges number) and guarantees continuity in successive u f 2 edges. By replicating this stochastic reconstruction process, a forest of n reconstructed ADN, denoted (d 1 , . . . , d n ), can be simulated. This forest represents uncertainties in ADN mapping.

  2.1.5. Neighbour state process: random pruning and branching During simulated annealing, neighbouring random states for digraph d are generated (see line 7 in algorithm 1). A neighbour state is based on random pruning (Fig. 1-b) and branching (Fig. 1-a) of the d graph as depicted on figure 1 on a virtual example.

Figure 1 :

 1 Figure 1: A neighbour state of tree d generated through random branching (a) or random pruning (b) processes. The two wide black segments are fixed parts of the network and the wide grey point corresponds to the tree root (outlet).

  denotes the cumulative length of f , d and f , respectively, within d edges. γ 1 , γ 2 are convergence tuning parameters that permit the differential adjustment of each of the objectives (e 1 , e 2 ) relative to the unique e max parameter.

  by 2 and normalised by the target catchment cumulative length. The α parameter was fixed such that exp ( and 0.1 values were fixed arbitrarily after trials and errors tests since automated fitted parameters for simulated annealing is still an open question[START_REF] Wright | Automating parameter choice for simulated annealing[END_REF]. There were also fixed according to suggestions from Aarts et al.(2003) for initial temperature related to the maximal cost between two neighbouring solutions and from[START_REF] Zhang | A simulated annealing with a new neighorhood structure based algorithm for high school timetabling problems[END_REF] in order to get a rejection rate of about 5-10 %.2.2. Evaluating the ADN simulations with drainage network similarity metricsThe previous subsection showed the way an ADN would be reconstructed given the following information: r, the known outlet of the network; l, the plot boundary lattice graph valued by elevation, and f , the set of fixed network reaches used to condition the reconstruction.

Figure 2 :

 2 Figure 2: ADN topographical similarity computation through empirical cumulative edge slope distribution difference (grey areas): example of small (a) and great (b) dissimilarities for ditree d1 and d2 compared to d ref

Figure 3 :

 3 Figure 3: The Muscat catchment data: a) the cultural plot lattice depicted on the 5-m hillshaded DTM; b) the actual drainage network in blue, with the reach width related to Strahler order ranging from 1 to 5; c) the plot lattice with the set f of known edges (blue: c BD-Topo for f2, red: remote sensing for f1)

Figure 4 :

 4 Figure 4: Simulated annealing algorithm 1 intermediary states: a) initial digraph and energy for k = 0; b) digraph and energy at iteration k = 30; c) digraph and energy at final iteration k = 309; d) final obtained ditree. On the top, the lattice edges are depicted in grey; the directed reconstructed graphs are indicated by black arrows. The fixed edges of f with a known direction ( c BD-Topo) are depicted in blue, and the fixed edges of f with an unknown direction (remote sensing) are depicted in red. On the bottom, the energy e is depicted as a function of k.

Figure 5 :

 5 Figure 5: Example of five reconstructed networks using a) r, b) f2, c) f1 +f2. All networks are depicted with a line width related to their Strahler orders.

Figure 6 :

 6 Figure 6: Similarity distribution of forests using a) r, b) f2, c) f1 + f2. On the top, the F θ (slope) curves are depicted for the actual network (red) and the reconstructed networks (black). On the bottom, the Strahler barycenter of the actual network is depicted by the black point from which red segments join the reconstructed network Strahler barycenters.

  , 2003. Local search in combinatorial optimization. Princeton: Princeton University Press, Ch. Simulated Annealing. Amini, J., Blais, M. S. J., Lucas, C., Azizi, A., 2002. Automatic road-side extraction from large scale imagemaps. International Journal of Applied Earth Observation and Geoinformation 4, 95-107.

  e max do

7: d new ← neighb(d) random neighbour state 8: e new ← energy(d new ) energy function 9: if e new < e then 10: d ← d new ; e ← e new 11:

Table 1 :

 1 Dissimilarities statistics

	length (%)	42.6	45.6	57.5
	Ld	23850.5 23868.5	24763
	ht	0.00101 0.00118 0.00394
	σ ht	0.00201 0.00217 0.00275
	hc	244.9	230.2	197.5
	σ hc	23.5	22.2	37.5

Table 2 :

 2 Welch two sample T test probability values

Author-produced version of the article published in International Journal of Applied EarthObservation and Geoinformation, 2011, 13(6), 853-862. The original publication is available at http://www.sciencedirect.com/ Doi : 10.1016/j.jag.2011.06.001

Author-produced version of the article published in International Journal of Applied Earth Observation and Geoinformation, 2011, 13(6), 853-862. The original publication is available at http://www.sciencedirect.com/ Doi : 10.1016/j. jag.2011.06.001 Author-produced version of the article published in International Journal of Applied Earth Observation and Geoinformation, 2011, 13(6), 853-862. The original publication is available at http://www.sciencedirect.com/ Doi : 10.1016/j.jag.2011.06.001 Author-produced version of the article published in International Journal of Applied Earth Observation and Geoinformation, 2011, 13(6), 853-862. The original publication is available at http://www.sciencedirect.com/ Doi : 10.1016/j.jag.2011.06.001

Author-produced version of the article published in International Journal of Applied Earth Observation and Geoinformation, 2011, 13(6), 853-862. The original publication is available at http://www.sciencedirect.com/ Doi : 10.1016/j.jag.2011.06.001


Yohan Thibault

Yukiko Kenmochi

Akihiro Sugimoto

Computing Upper and Lower Bounds of Rotation Angles From Digital Images

. A sequence of hinge angles is a set of particular angles determined by a digital image in the sense that any angle between two consecutive hinge angles gives the identical rotation of the digital image. We propose a method for obtaining the lower and the upper bounds of admissible rotation angles using hinge angles from a given Euclidean angle or from a pair of corresponding digital images.

Introduction

Rotations in the discrete plane are required in many applications for image computation such as image matching or construction of mosaic images [START_REF] Hansen | Real-time scene stabilization and mosaic construction[END_REF]. For the moment, the most popular method to estimate the rotation angle is to approximate the rotation matrix by minimizing errors [START_REF] Hansen | Real-time scene stabilization and mosaic construction[END_REF].

In the continuous plane, the Euclidean rotation is well-defined and possesses the property of bijectivity. This implies that for two angles γ 1 , γ 2 and a set of points A, if the Euclidean rotation with angle γ 1 applied to A gives the Email addresses: thibauly@esiee.fr (Yohan Thibault), y.kenmochi@esiee.fr (Yukiko Kenmochi), sugimoto@nii.ac.jp (Akihiro Sugimoto).

same result as the Euclidean rotation with angle γ 2 applied to A, then we have γ 1 = γ 2 .

In the discrete plane, however, two different points can arrive at the same grid point after a discretization of the Euclidean rotation (DER). Because of this reason, two different angles can give the same rotation result for a set A of grid points 1 . In other words, we can define a set of admissible rotation angles S such that any angle in S gives the same rotation result for A. Note that S depends on A. The two most interesting angles in S are the lower and the upper bounds because with only these two angles we can deduce the other angles in S. This paper aims to find these two angles from a given rotation angle or from two given corresponding sets of grid points. In order to identify the exact bounds, we should not involve any computation error. Thus, we work with discrete geometry tools which guarantee to avoid computation with real numbers. Moreover, because we assume that our data are discretized from continuous images of an object, we enforce the property that the discrete rotation2 between two different sets of grid points gives the same result as DER.

Some work on discrete rotations already exists. The most widely used discrete rotation in the beginning is the CORDIC algorithm [START_REF] Volder | The CORDIC trigonometric computing technique[END_REF]. The rotation angle is estimated with addition or subtraction using pre-computed values to achieve the required precision. It gives almost the same result as DER but an approximation of the angle. Andres described in [START_REF] Andres | Cercles discrets et rotations discr tes[END_REF], [START_REF] Andres | The quasi-shear rotation[END_REF] some discrete rotations such as the rotation by discrete circles, the rotation by Pythagorean lines or the quasi-shear rotation. Computation executed during these rotations are exact. However, because they preserves the bijectivity, they cannot give the same results as DER.

On the other hand, Nouvel and Rémila proposed in [START_REF] Nouvel | Incremental and transitive discretized rotations[END_REF] another discrete rotation based on hinge angles, which gives the same results as DER. It is known that a sequence of hinge angles is a set of particular angles determined by a digital image in the sense that any angle between two consecutive hinge angles gives the identical rotated digital image. This means that hinge angles correspond to the discontinuity of DER. Nouvel and Rémila showed that each hinge angle is represented by an integer triplet, so that any discrete rotation of a digital image is realized only with integer calculation. Because their algorithm gives the same rotation results as DER, we see that hinge angles represented by integer triplets give sufficient information for executing any digital image rotation.

In this paper, we propose a discrete method for finding the lower and the upper bounds of admissible rotation angles. Our method uses hinge angles because we can obtain the same result as DER and because they allow exact computations. The input data of our method is two sets of grid points where point correspondences across the two sets are known. The output is two hinge angles that give the lower and the upper bounds of the admissible rotation angles for the two sets. Note that a part of this work was presented in [START_REF] Thibault | Computing admissible rotation angles from rotated digital images[END_REF].

Discrete Rotation

It can appear strange that the Euclidean rotation used for the common task in geometric computation is problematic for many applications. Data are usually represented in the computer by integers or rational numbers. But the Euclidean rotation which requires sine and cosine functions is designed for real numbers. Therefore, the computation results given by a Euclidean rotation are, in most cases, represented by floating numbers which are approximations of real numbers. When an algorithm uses the Euclidean rotation for integer or rational data and then converts the obtained floating values into integer or rational numbers, precisions of these results may decrease. Another problem also arises for rotations in discrete space. It is well known that the Euclidean rotation is bijective and transitive. But when we convert the result obtained by the Euclidean rotation into a discrete space, it is easy to see that these two properties are lost [START_REF] Hersch | Raster rotation of bilevel bitmap images[END_REF][START_REF] Paeth | A fast algorithm for general raster rotation[END_REF]. The loss of these two properties leads to research on the discrete rotation.

There are two ways to compute a discrete rotation: using floating numbers and using only integers. The first way, in most cases, is easiest and allows us to use floating computation followed by the rounding function to obtain the set of grid points as the output. The main problem with the rounding function is the approximation due to the loss of the value after the decimal point. This approximation leads to lack of precision during computation. The second way does not have this problem, but avoiding floating numbers implies that sine and cosine functions should not be used. Computing rotations without trigonometrical functions requires development of a new method, which is a tough problem.

Hinge Angles

We consider grid points in 2 as the centers of pixels and rotate them in such a way that the rotation center has integer coordinates. Hinge angles are particular angles that make some points in 2 rotated to points on the frontier between adjacent pixels. In this section, we give the definition of hinge angles and their properties related to Pythagorean angles.

Definition of Hinge Angles

Let x = (x, y) be a point in Ê 2 . We say that x has a semi-integer coordinate if x + 1 2 ∈ or y + 1 2 ∈ . The set of points each of which has a semi-integer coordinate is called the half-grid and is denoted by H . Thus, H represents the set of points on the frontiers of all pixels whose centroids are points in 2 .

Definition. 1 An angle α is called a hinge angle if at least one point in 2 exists such that its image by the Euclidean rotation with α belongs to H .

Because H can be seen as the discontinuity of the rounding function, hinge angles can be regarded as the discontinuity of the discretized Euclidean rotation. More simply, hinge angles determine a transit of a grid point from a pixel to its adjacent pixel during the rotation.

The following theorem is important because it shows that we can represent every hinge angle with three integers.

Theorem. 2 ([7]

) An angle α is a hinge angle for a grid point (P, Q) ∈ 2 if and only if there exists K ∈ such that

2Q cos α + 2P sin α = 2K + 1. (1)
Geometrically, a hinge angle α is formed by two rays going through (P, Q) and a half-grid point (K + 1 2 , λ) where the two rays share the origin as their endpoints as shown in Figure 1 (left). This theorem indicates that all calculations related to hinge angles can be done only with integers. Hereafter, α indicates a hinge angle. We denote by α(P, Q, K) the hinge angle generated by an integer triplet (P, Q, K).

Setting λ = P 2 + Q 2 -(K + 1
2 ) 2 , we easily derive the following equations from (1) and Figure 1 (left),

cos α = P λ + Q(K + 1 2 ) P 2 + Q 2 , sin α = P (K + 1 2 ) -Qλ P 2 + Q 2 . ( 2
)
Fig. 1. A hinge angle α(P, Q, K) (left) and four symmetrical hinge angles (right).

Note that we have a case where a half-grid point is (λ, K + 1 2 ) instead of (K + 1 2 , λ). In such a case, the above equations become

cos α = Qλ + P (K + 1 2 ) P 2 + Q 2 , sin α = P λ -Q(K + 1 2 ) P 2 + Q 2 . ( 3
)
The symmetries on hinge angles are important because it allows us to restrict rotations in the first quadrant of the circle such that α ∈ [0, π 2 ].

Corollary. 3 Each triplet (P, Q, K) corresponds to four symmetrical hinge angles α + πk 2 where k = 0, 1, 2, 3.

Figure 1(right) gives an example of Corollary 3. In order to distinguish the case of (K + 1 2 , λ) from that of (λ, K + 1 2 ), we change the sign of K; we use α(P, Q, K) for the case of (K+ 1 2 , λ), and α(P, Q, -K) for the case of (λ, K+ 1 2 ). Because the symmetries allow us to restrict α to the range [0, π 2 ], as mentioned above, we may assume that K is positive.

Properties Related to Pythagorean Angle

Because hinge angles are strongly related to Pythagorean angles, certain properties of Pythagorean angles are required to prove some properties of hinge angles. Thus, we first give the definition of Pythagorean angles and their properties.

Definition. 4 An angle θ is called Pythagorean if and only if both its cosine and sine belong to the set of rational numbers É.

We can deduce from Definition 4 that a Pythagorean angle θ is represented by an integer triplet (a, b, c) such that

cos θ = a c , sin θ = b c . (4)
In the following, θ indicates a Pythagorean angle. The lemma below for Pythagorean angles is well known. Theorem. 6 Let E h be the set of hinge angles and E p be the set of Pythagorean angles. Then we have E h E p = ∅.

Proof. Assume that there exists an angle α such that α ∈ E h and α ∈ E p . Since α ∈ E p , we can find an integer triplet (a, b, c) generating α where gcd(a, b, c) = 1. By substitution of (4) in (2), we obtain

2 Qa + P b c = 2K + 1, (5)
from which we derive 2 Qa+P b c ∈ . Because we know that c is odd according to Lemma 5, we obtain Qa+P b c ∈ . However, this contradicts the fact that for any pair n, m ∈ , we never have 2n = 2m + 1. Therefore α cannot belong to E h and E p simultaneously.

This theorem shows that it is not possible to rotate a point (i, j) ∈ 2 to a point (x, y) such as x = i + 1 2 , y = j + 1 2 , if the angle of the rotation is a hinge angle.

The next theorem shows an interesting relation between hinge angles and Pythagorean angles.

Theorem. 7 ( [START_REF] Nouvel | Rotations discr tes et automates cellulaires[END_REF]) Let θ be a Pythagorean angle and α be a hinge angle. The angle α ′ = α + θ is a hinge angle.

Computing the Lower Bound Rotation Angle from a Pythagorean Angle

Because we work here in the discrete space, the rotation of a grid point by two different angles can give the same result. Namely, two different angles give the same result after the rotation of a grid point followed by discretization. Generally there exists a range of angles in which the same result is obtained. We thus define admissible rotation angles, abbreviated hereafter by ARA, to represent this range of angles.

In this section, we propose a method for computing the lower bound α inf of ARA for a given digital image from a given angle. Note that with minor modifications, this method can also find the upper bound α sup of ARA. Thus applying any rotation to the given digital image with an angle between α inf and α sup gives the same result. We note that both α inf and α sup are hinge angles.

Our input is a Euclidean angle. However, we can replace it by a Pythagorean angle because there exists a method with linear time complexity O(m) to approximate a given Euclidean angle with a Pythagorean angle with a precision of 1 10 m [START_REF] Anglin | Using Pythagorean triangles to approximate angles[END_REF], where m is a fixed integer that represents the quality of approximation. Below, we assume that we are given a Pythagorean angle, as the one in [START_REF] Nouvel | Incremental and transitive discretized rotations[END_REF].

Nouvel and Rémila presented a method for computing all possible hinge angles for a grid point or a pixel in a digital image [START_REF] Nouvel | Incremental and transitive discretized rotations[END_REF]. Their method can be used for finding our interesting hinge angle which is the lower bound of the admissible rotation angles. Its time complexity is O(n log(n)) where n is the number of all hinge angles for a given grid point. Note that n depends on the coordinates of the grid point.

Computing the Lower Bound Rotation Angle for a Grid Point

For each grid point p = (P, Q) ∈ 2 , there are less than n = ⌊ √ P 2 + Q 2 + 1 2 ⌋ different hinge angles in each quadrant [START_REF] Nouvel | Incremental and transitive discretized rotations[END_REF]. We can compare in magnitude any pair of hinge angles. This means that we have a totally ordered set {α(P, Q, K 1 ), α(P, Q, K 2 ), ..., α(P, Q, K n )} of hinge angles in the ascending order where K i ∈ . Given a Pythagorean angle θ, in order to find the lower bound rotation angle α(P, Q, K i ) such that α(P, Q, K i ) < θ < α(P, Q, K i+1 ), we use a tree structure. The binary search allows us to find α(P, Q, K i ) in O(log(n)), providing that we can compare a hinge angle with a Pythagorean angle in a constant time. The algorithm is described in Figure 2.

Function: Find the lower bound rotation angle for a point Input (Point p(P, Q), Pythagorean angle θ) The following theorem shows that the comparison between a hinge angle and a Pythagorean angle is executed in a constant time.

Output (α(P, Q, K)) var K max = ⌊ √ P 2 + Q 2 -1⌋; var K min = 0; var K = ⌊ Kmax+K min 2 ⌋; While (K max -K min = 1) if (α(P, Q, K) > θ) K max = K; else K min = K; K = ⌊ Kmax+K min 2 ⌋; end while return α(P, Q, K);
Theorem. 8 Let α be a hinge angle and θ be a Pythagorean angle. We can check whether α > θ in a constant time with integer calculation.

Proof. Let α(P, Q, K) be a hinge angle in [0, π 2 ] and θ(a, b, c) be a Pythagorean angle in [0, π 2 ]. From ( 2) and ( 4), we obtain

cos α -cos θ = Q(K + 1 2 ) + P λ P 2 + Q 2 - a c .
If θ is greater than α, cos αcos θ > 0. Thus

cQ(2K + 1) -2a(P 2 + Q 2 ) > -2cP λ. (6)
Since we know that c, P, λ are positive, the right-hand side of ( 6) is always negative. Thus, if the left-hand side of ( 6) is not negative, then θ > α. Otherwise, we take the square of ( 6), so that we only have to check whether the following inequality holds:

cQ(2K + 1) -2a(P 2 + Q 2 ) 2 < 4c 2 P 2 λ 2 . ( 7
)
Note that because

λ = P 2 + Q 2 -(K + 1
2 ) 2 , we see that 4λ 2 in the righthand side of (7) contains only integer values. Therefore, we can verify [START_REF] Nouvel | Incremental and transitive discretized rotations[END_REF] with integer calculation. If it is true, θ > α; otherwise, α > θ. Note that because of Theorem 6, it is impossible to obtain θ = α.

We claim that this comparison of a hinge angle and a Pythagorean angle is executed in a constant time because even in the worst case, we only have to check two equations ( 6) and [START_REF] Nouvel | Incremental and transitive discretized rotations[END_REF].

We mention the importance of the rotation with angle π 2 and its multiplications. In fact, if the angle of a rotation is equal to π 2 , π, 3π 2 , we just have to flip x and/or y-coordinates by changing their signs. It gives the justification that we can restrict the input angle θ to 0 < θ < π 2 .

Computing the Lower Bound Rotation Angle for a Set of grid Points

In this subsection, we present an algorithm for computing the lower bound rotation angle from a given Pythagorean angle θ for a digital image consisting of m grid points A. The output is a triplet of integers that represents the lower bound rotation angle for A. We note that the lower bound rotation is a hinge angle.

The algorithm computes all hinge angles for all points in A, and sorts them to keep the largest one. More precisely, we first compute the lower bound rotation angle for the first point of A, and store it as a reference. Then, we compute the lower bound rotation angle for the second point in A and compare it with the reference to keep the larger one. After repeating this procedure for all points in A, our algorithm returns the lower bound rotation angle α such that α < θ.

The time complexity of this algorithm is O(m log(n)) because we call m times the binary search (Figure 2) whose time complexity is O(log(n)). Figure 3 illustrates our algorithm. As shown in Theorem 9, the comparison between two hinge angles is realized in a constant time, so that our algorithm does not change the global complexity.

Theorem. 9 Let α 1 , α 2 be two hinge angles. We can check whether α 1 > α 2 in a constant time and with integer calculation.

Proof. Let α 1 (p, q, k) and α 2 (r, s, l) be two hinge angles in [0, π 2 ]. From (3) we obtain

cos α 1 -cos α 2 = p(k + 1 2 ) + qλ 1 p 2 + q 2 - r(l + 1 2 ) + sλ 2 r 2 + s 2 . ( 8
)
If α 2 is greater than α 1 , cos α 1 -cos α 2 > 0. Thus (r 2 + s 2 )p(2k + 1) -(p 2 + q 2 )r(2l + 1) > 2(p 2 + q 2 )sλ 2 -2(r 2 + s 2 )qλ 1 . (9)
If the left-hand side of ( 9) is negative and the right-hand side of ( 9) is positive, then α 1 > α 2 . If the left-hand side of ( 9) is positive and the right-hand side of ( 9) is negative, then α 2 > α 1 . We can easily check the signs of the lef-hand side and the right-hand side of ( 9) with integer computation. Note that p, q, k, r, s, l are all positive, and that (2(p 2 + q 2 )sλ 2 ) 2 and (2(r 2 + s 2 )qλ 1 ) 2 contain only integer values.

If the signs of the left-hand side and the right-hand side of ( 9) are the same, we first compute the square of each side and then compare the values to identify which is the greater. For simplicity, we assume that the signs of the both sides of ( 9) are positive, and let A = (r 2 + s 2 )p(2k + 1), B = (p 2 + q 2 )r(2l + 1), C = (r 2 + s 2 )q and D = (p 2 + q 2 )s. Now ( 9) is rewritten by

A -B > 2Dλ 2 -2Cλ 1 . (10)
Then we take the square of equation ( 10) to obtain

(A -B) 2 -4(C 2 λ 2 1 + D 2 λ 2 2 ) > -8CDλ 1 λ 2 . (11)
If the sign of the left-hand side of ( 11) is positive, we can deduce that α 2 > α 1 .

Otherwise, taking the square of each side gives us

(A -B) 2 -4(C 2 λ 2 1 + D 2 λ 2 2 ) 2 < 64C 2 D 2 λ 2 1 λ 2 2 . ( 12
)
We note that we can easily verify whether (12) is satisfied with integer computation alone. If (12) is true, α 1 < α 2 ; otherwise α 2 < α 1 . The same logic can be applied to the case where the signs of the both sides of (9) are negative.

This comparison of a pair of hinge angles is executed in a constant time because in the worst case, we only have to check three equations ( 9), ( 11) and (12).

Digital Image Rotation by a Lower Bound Rotation Angle

This section uses the results obtained in Section 4 to present an algorithm for rotating a digital image with a given lower bound rotation angle.

It is already proved in [START_REF] Nouvel | Incremental and transitive discretized rotations[END_REF] that we can obtain the same result as the DER with respect to the original rotation angle. Note that our input is a lower bound rotation angle and the input of the algorithm presented in Figure 2 is a Pythagorean angle. In spite of this difference, we can apply the same algorithm thanks to Theorem 9, since we are looking for K, which gives the arriving pixel (K, ⌊λ + 1 2 ⌋), for each pixel (P, Q). The algorithm is presented in Figure 4. It supposes that the center of rotation is the origin.

For each point, our algorithm calls the binary search (Figure 2) to find the corresponding hinge angle, which designates its new position. If we consider n as the biggest coordinate of all points in A, we can assume that there are less than 4n 2 points in A. Thus we can conclude that the complexity of our algorithm is O(n 2 log(n)). The first advantage of our method is that it does not require any floating number calculation. The second advantage is that the exact rotation of the digital image is obtained with only an integer triplet. We need neither matrices nor angles for realizing the rotation.

Obtaining Admissible Rotation Angles from Two Digital Images

Let us assume that a set of grid points in the first image and its corresponding set in the second image are given: A = {p 1 , p 2 , ..., p l } and B = {q 1 , q 2 , ..., q l } are given where p i corresponds to q i . Given A and B, we obtain a hinge angle pair {α inf , α sup }. This pair of hinge angles is the lower and the upper bounds of the ARA. Therefore each γ such that α inf ≤ γ < α sup is consistent with the point correspondences between A and B. Hereafter, we assume that A is the original point set and B is the rotated point set by angle γ. In this section, we show how to obtain the ARA from A and B.

Setting Rotation Centers

For any rotation, we need to set a rotation center. Without loss of generality, we may choose any grid point in a digital image for the rotation center. Assuming that rotation centers for A and B are p 1 and q 1 respectively, we define two translation functions T A and T B such that

T A (p i ) = p i -p 1 , T B (q i ) = q i -q 1 ,
for all p i ∈ A, q i ∈ B. We can regard the origin as the rotation centers after these translations. Hereafter, we assume that these translations have been already applied in order to obtain A = {p 1 , p 2 , ..., p l } and B = {q 1 , q 2 , ..., q l }.

Computing Lower and Upper Bounds of Rotation Angles from Two Corresponding Point Pairs

In this subsection, we consider a special case of two corresponding point pairs.

We let A = {p 1 , p 2 } and B = {q 1 , q 2 } where p i = (P i , Q i ) and q i = (R i , S i ).

We then define a circle C (p 2 ) going through p 2 whose center is p 1 . Thus the radius of C (p 2 ) is r = d(p 1 , p 2 ) where d(p 1 , p 2 ) is the Euclidean distance between p 1 and p 2 .

Let us define the half-grid H (q 2 ) around q 2 : Fig. 5. The corners of H (q), namely, four corners of a pixel around q.

H (q 2 ) = {(x, y) ∈ H :

S 2 - 1 2 ≤ y ≤ S 2 + 1 2 if x = R 2 + 1 2 , R 2 - 1 2 ≤ x ≤ R 2 + 1 2 if y = S 2 + 1 2 }.
We set p 1 and q 1 to be the rotation centers. Then, we need to detect intersections between C (p 2 ) and H (q 2 ) in order to find a hinge angle pair. We thus investigate which corners of H (q 2 ) are inside of C (p 2 ).

Setting four corners of H (q 2 ) to be

C 1 (q 2 ) = (R 2 -1 2 , S 2 -1 2 ), C 2 (q 2 ) = (R 2 -1 2 , S 2 + 1 2 ), C 3 (q 2 ) = (R 2 + 1 2 , S 2 + 1 2 ), C 4 (q 2 ) = (R 2 + 1 2 , S 2 - 1
2 ) as shown in Figure 5, we define a binary function F :

F (C i (q 2 )) =    1 if C i (q 2 ) is inside of C (p 2 ), 0 otherwise.
In order to obtain F (C i (q 2 )) with integer calculation, we compare each of (2(R 2 ± 1 2 ), 2(S 2 ± 1 2 )) 2 with (2r) 2 . Note that we may assume that C (p 2 ) and H (q 2 ) always intersect with each other. This is because no intersection between C (p 2 ) and H (q 2 ) indicates that p 2 and q 2 are not corresponding.

The following lemmas are needed to prove Theorem 12 below.

Lemma. 10 For a circle C (p 2 ) centered on p 1 , any

C i (q 2 ) = (R 2 ± 1 2 , S 2 ± 1 2 ) cannot be on C (p 2 ) for i ∈ {1, 2, 3, 4} where R 2 , S 2 ∈ .
Proof. Let r be the radius of C (p 2 ). Because C (p 2 ) goes through p 2 , r 2 ∈ . Let us assume that C i (q 2 ) is on C (p 2 ). Thus the Euclidean distance between the origin and C i (q 2 ) is equal to r. This indicates that r 2 = (R 2 ± 1 2 ) 2 +(S 2 ± 1 2 ) 2 . However, this contradicts the above fact that r 2 ∈ .

Lemma. 11 Let D be a line that belongs to H . If C (p 2 ) is a circle centered on p 1 . Then, the number of distinct intersections of C (p 2 ) and D is two or zero.

Proof. Let p 2 = (P 2 , Q 2 ) ∈ and the equation representing D be x = i + 1 2 (i ∈ ). Letting (x 1 , y 1 ) be the coordinates of the intersecting point of D and C (p 2 ). Then we have

   x 2 1 + y 2 1 = P 2 1 + Q 2 1 , x 1 = i + 1 2 .
From these equations, we obtain

y 2 1 = P 2 1 + Q 2 1 -(i + 1 2 ) 2 . Because P 2 1 + Q 2 1 - (i + 1
2 ) 2 does not belong to , y 1 cannot be equal to 0. Therefore there is two distinct solutions for y

1 if P 2 1 +Q 2 1 -(i+ 1
2 ) 2 > 0; no solution otherwise. Similar discussion can be done for the case of y = j + 1 2 .

Theorem. 12 If two points p 2 and q 2 are corresponding, the circle C (p 2 ) and the half-grid H (q 2 ) always have two or four distinct intersections.

Proof. In general, if a circle intersects with a square and the center of the circle is not inside of the square, we have 1, 2 or 4 intersections. Having just one intersection means that the circle goes through one corner of the square or the circle is tangential to a half-grid. Lemma 10 shows that the circle cannot goes through one corner. Lemma 11 shows that the circle cannot be tangential to any half-grid. Therefore we have only 2 or 4 intersections.

From Theorem 12, we always have two or four distinct intersections between C (p 2 ) and H (q 2 ), and we see that there are five cases corresponding to different possibilities to have 0,1,2 or 3 corners inside of C (p 2 ), as illustrated in Figure 6.

Remark. 13

The cases represented in (a) and (c) in Figure 6 never happen unless two conditions are satisfied.

(1) q 2 is on the x-axis or on the y-axis.

(2) The radius r of C (p 2 ) is sufficiently close to a half-integer.

The case of Figure 6(c) is problematic because it generates two distinctive ranges of admissible rotation angles. These two ranges are symmetric with respect to the x-axis or the y-axis. However, this case rarely occurs as explained below.

Supposing that q 2 is neither on the x-axis nor the y-axis, we see R 2 or S 2 is not zero. We assume that they are positive. In the first quadrant, the ycoordinate (respectively x-coordinate) of points in C (p 2 ) is strictly decreasing with respect to x (respectively y). Thus it cannot intersect twice with a line parallel to the x-axis (respectively y-axis).

Supposing that q 2 is on the y-axis (respectively on the x-axis), we see that the distance between the origin and C 1 (q 2 ) (respectively C 3 (q 2 )) is greater than r.

Thus we have r < (R 2 -1 2 ) 2 + ( 1 2 ) 2 (respectively r > (R 2 + 1 2 ) 2 + ( 1 2 ) 2 ). (a) ΣF (C i (q 2 )) = 0 (b) ΣF (C i (q 2 )) = 1 (c) ΣF (C i (q 2 )) = 2 (d) ΣF (C i (q 2 )) = 2 (e) ΣF (C i (q 2 )) = 3
Fig. 6. Illustration of cases ΣF (C i (q 2 )) = 0,1,2 or 3.

Letting ǫ = r -(R 2 -1 2 ) (respectively ǫ = r -(R 2 + 1 2
)) where ǫ < 1 2 , we obtain ǫ(8r -4ǫ) < 1 (in both cases). We experimentally observe that as far as r ≤ 10 5 , ǫ(8r -4ǫ) < 1 never holds. On the other hand, if r is sufficiently large enough to have 8r -4ǫ ≈ 8r then we have ǫ < 1 8r , which means that r is sufficiently close to a half-integer. We can thus conclude that (2) in Remark 13 is satisfied only if r is sufficiently large.

As explained above, the case of Figure 6(c) occurs extremely rarely. Thus, we take no account of this case from now on, which causes no problem from the practical point of view.

The main function of our algorithm for finding the lower and the upper bounds of admissible rotation angles, consists of three steps. The first step sets the rotation center at p 1 and q 1 , as described in Section 6.1. The second step computes which corners are inside of C (q 2 ) and then compute the index I q 2 = i 2 i × F (C i (q 2 )). Therefore we can easily identify which corners are inside of C (p 2 ) from I q 2 . The third step calls a function that returns hinge angles corresponding to I q 2 . There exist fourteen possible values for I q 2 from 0 till 15 except for 5 and 10. Note that geometrically I q 2 can be neither 5 nor 10. The value 15 of I q 2 implies an error such that all corners are inside of C (q 2 ). Since Table 1 Corners of H (q 2 ) inside of C (p 2 ) and ARA.

Value of I q

2 α inf α sup I q 2 = 1 or I q 2 = 14 α(P 1 , Q 1 , R 2 -1) α(P 1 , Q 1 , 1 -S 2 ) I q 2 = 2 or I q 2 = 13 α(P 1 , Q 1 , R 2 -1) α(P 1 , Q 1 , -S 2 ) I q 2 = 3 or I q 2 = 12 α(P 1 , Q 1 , 1 -S 2 ) α(P 1 , Q 1 , -S 2 ) I q 2 = 4 or I q 2 = 11 α(P 1 , Q 1 , R 2 ) α(P 1 , Q 1 , -S 2 ) I q 2 = 6 or I q 2 = 9 α(P 1 , Q 1 , R 2 -1) α(P 1 , Q 1 , R 2 ) I q 2 = 7 or I q 2 = 8 α(P 1 , Q 1 , R 2 ) α(P 1 , Q 1 , 1 -S 2 ) I q 2 = 0 and R 2 = 0 α(P 1 , Q 1 , 1 -S 2 ) α(P 1 , Q 1 , 1 -S 2 ) I q 2 = 0 and S 2 = 0 α(P 1 , Q 1 , R 2 -1) α(P 1 , Q 1 , R 2 -1)
I q 2 whose value is 0 corresponds to the case ΣF (C i (q 2 )) = 0, we should verify whether H (q 2 ) really intersects with C (p 2 ). Note that for the other values for I q 2 , we can make a pair (d, e) such that d + e = 15. The two indices of each pair design the same pair of lower and upper bounds of ARA. Table 6.2 gives the corresponding lower and upper bound rotation angles for each value of I q 2 . Each step of this algorithm has the constant time complexity. Thus the global complexity of this algorithm is also O(1).

Incremental Computing Lower and Upper Bounds of Rotation Angles

In general, the corresponding point sets contain more than two points. Therefore, in this section, we extend our algorithm in Section 6.2 to two sets of corresponding point pairs, A and B, each of which has l points where l > 2.

To simplify the notation, we denote by ARA(p i , q i ) = (α i inf , α i sup ) the pair of angles that gives the lower and the upper bounds of admissible rotation angles for the pair of points (p i , q i ). Note that α i inf , α i sup are hinge angles. ARA(A n , B n ) denotes the two most restrictive angles for all points i such as i ≤ n. We recursively define it by

ARA(A n , B n ) = ARA(A n-1 , B n-1 ) ARA(p n , q n ).
A new algorithm handles all points incrementally. This algorithm is divided into two parts. The first part is to initialize the algorithm by computing ARA(p 2 , q 2 ). The second part computes ARA(A i , B i ) for i = l. Note that ARA(p 1 , q 1 ) cannot be computed because p 1 and q 1 are the centers of the rotation.

The time complexity of this algorithm is O(l). As explained in Section 6.2, the function giving a pair of lower and upper bound rotation angles from a pair of points is realized in a constant time O(1). Moreover, as explained in Section 3, we can compare two hinge angles in a constant time O(1). Therefore, the computation of this algorithm for l points takes the time complexity of l × (O(1) + O(1)) = O(l) as a whole. Figure 7 gives an example of the incremental algorithm for two sets of three points. Given input data of the algorithm as shown in Figure 7 (A), we first obtain the result of the translation described in Section 6.1 as illustrated in (B). We then compare, for each pair of points (p i , q i ) with i ≥ 2, the distance of p i from the origin with that of each corner from H (q i ) to deduce the corresponding hinge angle as explained in Section 6.2. Finally, we obtain (D) which shows the intersection of all ARA(p i , q i ) obtained in (C).

Experimental Evaluation using Synthetic Data

We evaluated our algorithm using synthetic data and verified the quality of obtained admissible rotation angles. To test our algorithm, we need two sets of points. We randomly generated the set F of 100 floating points in a 200 × 200 square. The first set I is obtained by discretizing F . The second set I ′ is obtained by applying the DER with angle θ = 50 • to the first set I (Figure 8). than 0.1 • while the range is reduced to 0.02 • after twenty pairs of points. Then the precision increases slowly. This shows that the precisions of admissible rotation angles acquired after twenty pairs of points are not significant in this particular case.

Figure 9 also indicates that most pairs of points give the admissible rotation angle range smaller than 1 • (for example, pair #39). Pair #34 and pair #66, however, give a range greater than 4 • . In fact, the maximum range for a pair of points is directly related to the distance between the center of rotation and the pair of points. We denote by d the distance between the rotation center and the pair of points. Then, the maximum difference between the lower and the upper bound rotation angles is sin -1 √ 2 d . For the pair #34, d = √ 173; thus the maximum range is approximately 6.1 • . For the pair #39, d = √ 16505 and thus the maximum range is approximately 0.63 • . This is consistent with our experimental result.

We cannot expect that just by applying our algorithm only to the farthest pairs of points from the rotation center, we will obtain the exact bounds of admissible rotation angles. This can be confirmed by the pair #34, for example. We see that the range is greater than 4 • but the upper bound hinge angle gives a close approximation to the rotation angle. On the other hand, by keeping only the farthest pairs of points, our algorithm might become faster. But we need an additional procedure for determining which pairs of points have to be kept. If we set distance δ and keep only the pairs of points whose distance from the rotation center is greater than δ, we have to compute the distances of all pairs of points from the rotation center. Then the complexity required by this approach is still O(n). Accordingly, we conclude that keeping only the farthest pairs of points is a does not contribute toward the complexity improvement.

Discussion on practical application to digital images

In the previous section, we start with grid points in the discrete plane. In other words, all coordinates of the points are integers. But in reality, points in the Euclidean space are not represented by integers but real numbers. Thus we start here with the set of floating points to see how our algorithm works. Then we test our algorithm with real data acquired by a digital camera.

Synthetic data

To test our algorithm using synthetic data more similar to real data, we applied the Euclidean rotation with angle θ = 50 • directly to the set F that is used in Section 7. Then we discretized the rotated image to obtain the second set I ′′ . To the two sets I and I ′′ , we applied our algorithm. We note that the set I comes from Section 7.

Figure 10 shows the obtained hinge angles. Because we started with the set of floating points, we observe some errors in bounding admissible rotation angles. In fact, we see that some pairs of points, the pairs #64 and #73 for instance, do not contain θ.

Figure 11 (left) explains how errors arise in rotating floating points. Let assume that the point q 1 is obtained after the Euclidean rotation of the floating point p 1 with angle γ. Let p 2 and q 2 denote the discretization of p 1 and q 1 . If we apply our algorithm to the pair of points {p 2 , q 2 }, we obtain ARA(p 2 , q 2 ) = (α inf , α sup ). In this case, γ does not belong to the interval [α inf , α sup ] (Figure 11 (right)). We see that this is caused by discretization of the floating points. We can give a bound for this error. In fact, this bound directly depends on the distance between the rotation center and the pair of points. If we denote by d this distance, the maximum error is equal to sin -1 √ 2 2d . Note that the pair of points {p 3 , q 3 } and the angle ρ bring the same problem.

To avoid this problem, we modify our algorithm for example as follows. We keep all the pairs of bound rotation angles determined by all pairs of points. After sorting all lower and upper bound rotation angles into two lists L inf , L sup , we remove all hinge angles α inf from L inf such that for ∃α sup ∈ L sup , α inf > α sup . We also remove all hinge angles α sup from L sup such that for ∃α inf ∈ L inf , α sup < α inf . With this procedure, we can guarantee that the lower bound hinge angle is smaller than the rotation angle applied to the points and that the upper bound hinge angle is larger than the rotation angle. The complexity O(n log n) is required in sorting the remaining hinge angle pairs, which does not increase the computational cost of the algorithm as a whole.

Another way to avoid the problem is to compute the smallest distance d of all the points from the rotation center and then compute a Pythagorean angle θ satisfying θ > sin -1 √ 2 2d . We then add (subtract) θ to the upper (lower) bound rotation angle obtained by our algorithm in Section 6. As a result we obtain the upper and the lower rotation angles that define the interval accurately including the true rotation angle.

Both methods return a valid ARA. It is preferable to use the second method when there are few pairs of points because it does not remove any hinge angles. The first method is preferable when data contains many pairs of points in correspondence because the obtained ARA is more restrictive than the ARA obtained by the second method.

Real data

We applied our algorithm to see its practical usefulness. We used a turntable that is rotated with respect to the vertical axis with respect to a digital image plane. The precision of rotations in control is 10 -3 degrees.

We put a toy block on the turntable and then took its image using a standard digital camera where the camera was fixed so that its optical axis is parallel with the rotation axis (Figure 12(a)). Next we rotated the turntable with the angle of 44.99 • and then took another image of the toy block by the fixed camera (Figure 12(b)).

We manually selected five points in the first image and their corresponding points in the second image. Then we applied our algorithm to the five pairs of corresponding points.

Our algorithm might return the empty result (Figure 12(c)). This is because no intersection between C (p i ) and H (q i ) was found for all i = 1, 2, . . . , 5.

In the case of real data, we cannot always guarantee to detect correct corresponding pairs of points even manually. This indicates that a problem different from the discretization problem (see Section 8.1) arises 3 . Namely the maximum difference of distances from the rotation center between two points in correspondence can become greater than √ 2. This is illustrated in Figure 13 where points p and q are in correspondence and ǫ d = |d pd q | is greater than √ 2 where d p , d q are respectively the distances from p, q to the origin. We see that no intersection exists between the circle C (p) going through p and H (q).

To avoid this problem on corresponding points, we can change the resolution of the images. Namely, we degrade the image resolution until we find an intersection between C (p) and H (q). Because of the change in image resolution, we cannot accept the obtained hinge angles as they are. Instead, we add (subtract) a Pythagorean angle θ from the obtained upper (lower) hinge angle, in order to, let assume that d p < d q and that H (q) is defined for the image whose resolution is degraded with 2 -n from the original one. Then we have to choose θ satisfying θ > n √ 2 2dp .

Figure 12 (d) shows the results obtained by the modification above to the same data. We can see that the real rotation angle (44.99 • ) is included in the ARA obtained by our algorithm with this modification. Accordingly, we can conclude that our modification is effective.

Conclusion

In this paper, we have introduced problems of rotations in the discrete space. Because rotations are mathematically defined for the continuous space, it is necessary to develop new operational tools for such discrete rotations.

Our problem was from two digital images to find all possible rotation angles to rotate the first image into the second. Namely, to find the admissible rotation angles between two digital images. Since we need to have only exact computation without any approximation. We decide to adopt hinge angles which allow us rotations using only integer computation. By using their properties, we have proposed a method for incrementally computing the lower/upper bound rotation angles from a Pythagorean angle for a set of discrete points. Based on this method, we have also proposed a linear algorithm for finding from two digital images the lower and the upper bound rotation angles representing the admissible rotation angles between them.

In addition, we have shown with experiments that our algorithm is efficient with synthetic data and gives consistent results. Then we also have done some experiments to real data. Results given by our algorithm were not consistent because of the loss of correct point correspondences in our input. We then explained why these problems appear and gave necessary modifications to obtain consistent results. These modifications did not increase the time complexity while keeping the integer computation.

The range size of ARA can be considered as the unreliability of the point correspondence. In other words, the larger the range of ARA between two points becomes, the less the reliability of their point corresepondence is. Therefore, it may be interesting if we can propose a new method for evaluating "good" point correspondences from a given pair of digital images by using our results.

Another future work is to extend our proposed method to the 3D space. Because in 3D hinge angles are not yet defined in 3D, the first direction we should take is to design 3D hinge angles and develop a 3D discrete rotation algorithm based on 3D hinge angles. Then we will improve such a rotation algorithm in order to find admissible rotation angles in 3D as we did for 2D in this paper.

Lemma. 5

 5 Let (a, b, c) be an integer triplet generating a Pythagorean angle where |c| = max{| a |, | b |, | c |}. If gcd(a, b, c) = 1, then c is odd. If gcd(a, b, c) = i, then gcd( a i , b i , c i ) = 1 and the triplet of integers ( a i , b i , c i ) generates the same Pythagorean angle as (a, b, c).

Fig. 2 .

 2 Fig. 2. Function for the lower bound rotation angle for a point.

Function:

  Find the lower bound rotation angle for a digital image Input (Digital image A, Pythagorean angle θ) Output (hinge angle) var HA, HA temps \* hinge angle *\; HA = Find the lower bound rotation angle for a point (first point p 1 of A, θ); for each p ∈ A\{p 1 } HA temps = Find the lower bound rotation angle for a point (p, θ); if (HA < HA temps ) HA = HA temps ; end for Return (HA);

Fig. 3 .Fig. 4 .

 34 Fig. 3. Function for finding the lower bound rotation angle for a digital image.

Fig. 7 .

 7 Fig. 7. Running of the incremental algorithm.

Fig. 10 .

 10 Fig. 10. Result of our algorithm applied on a sets of floating points.

Fig. 11 .

 11 Fig. 11. Examples of errors for computing lower and upper bounds of rotation angles introduced by rotation of floating points.

Fig. 12 .

 12 Fig. 12. Result of our algorithm applied on real data.

Accordingly, DER is not bijective.

A discrete rotation is a rotation designed for the discrete space. It transforms a set of grid points into another set of grid points.

We assume here that we are to find ARA from given correspondences

Acknowledgment

This work was supported by grants from Région Ile-de-France.

We assumed that the point correspondences across the two sets I and I ′ are known.

Figure 9 shows the hinge angles obtained by our algorithm. The green and the red curves give the lower and the upper bound rotation angles for each pair of point in correspondence. We can see that lower bound rotation angles are always lower than θ and that upper bound rotation angles are greater than θ. The blue and purple lines give the lower and the upper bound rotation angles for points with the lower IDs than the point of interest. As we can see, after only ten pairs of points, the range of admissible rotation angles becomes less