Guillaume Brigot

Elise Colin-Koeniguer

Aurélien Plyer

Fabrice Janez

Adaptation and Evaluation of an optical flow method applied to co-registration of forest remote sensing images

The coregistration of heterogeneous geospatial images is useful in various remote sensing applications. Since the number of available data increases and the resolution improves, it is interesting to have an approach as automated, fast, robust and accurate as possible. In this paper, we present a solution based on optical-flow computation. This algorithm called GeFolki allows the registration of images in a non-parametric and dense way. GeFolki is based on a local method of optical flow derived from the Lucas-Kanade algorithm, with a multi-scale implementation, and a specific filtering including rank filtering, rolling guidance filtering and local contrast inversion. The efficiency of our coregistration chain is shown on radar, LIDAR and optical images on Remningstorp forest in Sweden. An analysis of the relevant parameters is investigated for several scenarios. Finally, we demonstrate the accuracy of our coregistration by proposing specific metrics for LIDAR/radar coregistration, and optics/radar coregistration.

I. INTRODUCTION

Coregistration of heterogeneous images is useful in various remote sensing image fusion applications, since one expects a gain from the synergy of sensors. Applications where data fusion is relevant are numerous, whether for land classification [START_REF] Wegner | Comparison of land cover classification using high-resolution terrasar-x and optical imagery[END_REF], for agriculture [START_REF] Qi | Synergy of optical and radar remote sensing in agricultural applications[END_REF], or forestry applications. In the case of the forest, many works illustrate the benefit of the combination of LIDAR and radar images [START_REF] Hyde | Exploring LiDAR RaDAR synergy predicting aboveground biomass in a southwestern ponderosa pine forest using LiDAR, SAR and InSAR[END_REF], [START_REF] Kaasalainen | Combining lidar and synthetic aperture radar data to estimate forest biomass: Status and prospects[END_REF], LIDAR and high resolution optical images [START_REF] Paris | A three-dimensional model-based approach to the estimation of the tree top height by fusing low-density lidar data and very high resolution optical images[END_REF], or radar and optical images [START_REF] Yitayew | Multisensor data fusion and feature extraction for forestry applications[END_REF], [START_REF] Hame | Improved mapping of tropical forests with optical and sar imagery, part i: Forest cover and accuracy assessment using multi-resolution data[END_REF], for the characterization of plant species or biomass assessment.

In all cases where the expected product is a map, geometrically aligning two or more images in order to combine pixels corresponding to the same objects is a crucial step of the fusion.

Most methods of remote sensing image coregistration are based either on geocoding, or on non-rigid image registration methods that use only the images as input.

In the case of geocoding, the accuracy of coregistration will be highly dependent on the availability and precision of both a DTM (Digital Terrain Model) and the orbit parameters [START_REF] Schubert | COSMO-SkyMed, TerraSAR-X, and RADARSAT-2 geolocation accuracy after compensation for earth-system effects[END_REF].

On the other hand, non rigid image registration without geocoding is widely investigated in various fields beyond the scope of remote sensing, for example in computer vision and medical imaging [START_REF] Jian | Robust nonrigid multimodal image registration using local frequency maps[END_REF]. In computer vision, video image co-registration has to meet the constraints of robustness and speed of execution, but often focuses on images taken from the same sensor with little delay in time. In medical imaging [START_REF] Gaens | Non-rigid multimodal image registration using mutual information[END_REF] or remote sensing, difficulty generally lies in the different nature of the images to compare. Moreover, the context of remote sensing is also changing today with larger quantities of time series data, and some time-sensitive applications require fast processing. This is the case for example for near real time change detection for rapid post disaster assessment, wildlife tracking, and surveillance across broad areas such as battlefields or border regions [START_REF] Coulter | Repeat station imaging for rapid airborne change detection[END_REF].

Most non-rigid registration methods are parametric methods, meaning that an assumption is made about a parametrized model that constrains the form of the expected deformations between processed images. Then a similarity function is optimized to find an approximation of a real underlying deformation [START_REF] Wong | ARRSI: automatic registration of remotesensing images[END_REF]. Among them, feature-based approaches establish a correspondence between a number of especially distinct points in images [START_REF] Li | Least squares image matching: A comparison of the performance of robust estimators[END_REF]. Selection of these points can use SIFT [START_REF] Goncalves | Automatic image registration through image segmentation and SIFT[END_REF], [START_REF] Fan | Registration of optical and SAR satellite images by exploring the spatial relationship of the improved SIFT[END_REF] or SURF methods [START_REF] Song | Remote sensing image registration based on retrofitted SURF algorithm and trajectories generated from lissajous figures[END_REF]. Other methods handle more complex features, such as segments, or use the shape descriptors [START_REF] Shi | Robust image registration using structure features[END_REF], [START_REF] Cao | A stable land cover patches method for automatic registration of multitemporal remote sensing images[END_REF].

The choice of an image similarity measure is a key point. One of the most widespread used for the registration of multimodality images is mutual information [START_REF] Chen | Mutual informationbased image registration for remote sensing data[END_REF]. Already used in remote sensing image coregistration, the main drawback of this measure is that it is quite time consuming. Instead of employing a similarity metric, [START_REF] Hu | Robust registration by rank minimization for multiangle hyper/multispectral remotely sensed imagery[END_REF] proposes to exploit a low rank constraint to jointly register multiple hyper-spectral images. Although such a model appears to be stable with respect to occlusions and imaging artifacts, it is not directly applicable in our case. In practice such a model requires a stack of multiple images and is further quite slow.

Another family of non-rigid coregistration methods are nonparametric. Among them, dense methods compute a displacement for every pixel in the image. They are particularly interesting in the case of very local deformation due for example to terrain elevation that has a lot of influence on high resolution images. Most of dense and non parametric methods belong to optical flow estimation. Optical flow is the pattern of apparent motion of objects in a visual scene caused by the relative motion between the sensor and the scene. Optical flow methods have been developed in a context where the constraints of speed and robustness to environmental effects have led to intensive efforts in producing algorithms that combine robustness, precision and high computing speed.

In this paper, we want to show how one optical flow method can undertake such a task of co-registration of heterogeneous images.

Optical flow has already been considered by the remote sensing community. [START_REF] Vogel | Optical flow for glacier motion estimation[END_REF] proposes a dense method to reconstruct the annual motion of glacier surfaces from ortho-photos.

In contrast, in this work we focus on the more general case of heterogeneous data, acquired by different sensor modalities, like the co-registration of LIDAR to SAR or optical images. In comparison to other existing methods developed for heterogeneous images often based on the use of mutual information (eg. [START_REF] Gaens | Non-rigid multimodal image registration using mutual information[END_REF]), our algorithm prevails clearly when the execution speed matters. Further, we demonstrate that our algorithm can be applied successfully to our needs, while providing added value compared to a method based on geocoding. We have already shown in [START_REF] Plyer | A New Coregistration Algorithm for Recent Applications on Urban SAR Images[END_REF] that an optical flow algorithm can successfully coregister SAR images in various conditions, including the demanding interferometric ones. In this paper we present a new version of this algorithm, called GeFolki, adapted to the coregistration of heterogeneous data, such as SAR-LIDAR images, or SAR-optical images. The acronym GeFolki stands for Geoscience Extended Flow Optical Lucas-Kanade Iterative.

A flow algorithm computes the displacement between two images. Most classical methods use partial derivatives of image intensity with respect to the spatial and temporal coordinates. They are based on two key points:

• First, an observation model where intensity remains the same from one image to another. This hypothesis is called the brightness constancy constraint. • Then, a linearization step in order to be able to solve the equation. Generally, it enables to write the variation of the associated intensity as a local gradient. In the case of heterogeneous images, brightness constancy is not respected because images have very different characteristics. It does not even make sense to compare SAR and photometric data with regard to intensity. Hence, the adaptation of such an optical flow algorithm to heterogeneous images is a challenging task. The solution proposed in this paper is to apply additional functions on the images to make them as similar as possible.

Moreover, the multi-scale strategy is used to estimate large displacements.

Our work is therefore based on the adaptation of a generic optical flow algorithm, that includes this multi-scale iterative strategy to compute the existing displacement between two remote sensing images of different types. The resulting new algorithm is called GeFolki. In section II, we present the previously mentioned GeFolki algorithm, detailing the different steps. Different results are presented in section III. Section IV assesses the underlying robustness of the result and deals with the parameter setting. In section V, we present a specific method to evaluate the performance of the results on the coregistration of images over forest lands, before concluding.

II. DESCRIPTION OF GEFOLKI

A. GeFolki position among optical flow algorithms

Let us consider the registration of two images I 1 and I 2 , defined on a 2D support S ∈ R 2 . The dense optical flow in computer vision is the pattern of apparent motion of objects. It corresponds to the displacement to find between both images. It is defined by u : x → u(x) ∈ R 2 . Those variables introduced are illustrated in Fig. 1. Most optical flow algorithms have a common starting point: it is to assume that pixel intensities remain constant from one frame to the next, what can be written:

I 1 (x) = I 2 (x + u(x)) (1)
This choice of observation model is called the brightness constancy. To compute the flow, the equality above needs usually a linearizing step. In most widely used techniques for optical flow computation, a change in intensity is written as:

I 2 (x + u(x)) = I 2 (x) + ∇I 2 • u (2)
Equation ( 2) is often called the linearized optical flow constraint. This constraint is based on the assumption that the displacements are small.

In the case of remote sensing images, almost all assumptions are violated: there is no temporal proximity between images, they are acquired under uncontrolled acquisition conditions and even by different sensors. The classical algorithms are thus not suited to model or compensate these effects. In this work we will show that with properly chosen modifications, the classical formulation can be adapted to these more general conditions.

As the optical flow system of Eq. 1 and 2 is underdetermined, optical-flow methods introduce in general additional conditions on the spatial distribution of the flow. Most prominently is to assume spatial smoothness in the flow field. This spatial additional condition can either be expressed globally or locally and relative methods can be classified respectively as local and global methods. Founding works are the global regularized framework of Horn and Schunck [START_REF] Horn | Determining optical flow[END_REF], and the local window-based approach of Lucas-Kanade [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF].

Local optical flow methods are often considered as more time and memory efficient, while global techniques are based on a richer prior formulation on the estimated motion flow [START_REF] Bruhn | Lucas/kanade meets horn/schunck: Combining local and global optic flow methods[END_REF]. Therefore, global methods can produce piecewise constant flow fields, as for example in the case of TV regularisation [START_REF] Zach | A duality based approach for realtime tv-l1 optical flow[END_REF]. However they have a higher computational cost due to global formulation of the problem. More recent approaches combine both concepts.

However in this article, our aim is not to produce the most efficient gradient-based optical flow algorithm, but to show that it is possible to adapt one of them for the coregistration of heterogeneous remote sensing images. Our choice was to adapt eFolki, a fast and robust local method, developed in the perspective of the trade-off between quality of the estimated flow and computational cost [START_REF] Plyer | Massively parallel Lucas Kanade optical flow for real-time video processing applications[END_REF]. EFolki is derived from the Lucas Kanade gradient-based method. Since eFolki has a remarkably simple and parallel structure, it is ideally suited to massively parallel computing architectures. Thus, the computation of the flow with our matlab implementation of GeFolki can be obtained in a few seconds for a 5,000 x 5,000 pixel image on a Intel(R) Core(TM) i7-2620M CPU @ 2.80GHz. A GPU implementation could provide this registration at the scale of hundreds of milliseconds.

Classical Lucas-Kanade algorithms assume that the flow is essentially constant in the local neighbourhood of the pixel under consideration, and solves the basic optical flow equations for all the pixels in that neighbourhood, by the least square criterion. Minimization is done by an iterative Gauss-Newton strategy based on first order Taylor expansion of the intensity of image around a previous guess of the optical flow u k , as illustrated in Fig. 2.

Figure 2. General minimization scheme by an iterative method. On the left: master image with a local window under consideration. In the middle: Two candidate regions in the slave image. The correct match is highlighted in red, another window with high similarity in blue. On the right: The matching cost is in general a non-convex function. Below: When started within sufficient proximity to the true minimum, gradient descent based optimization can recover the correct solution.

Modern Lucas-Kanade algorithms are not only iterative, but also multi-resolution. A pyramid of images is used to compute optical flow at varying scales following a coarse-to-fine strategy. This strategy enables to estimate large displacements despite the local constraint.

Our fast and robust algorithm can also be categorized as being of Lucas Kanade type. Thus, it belongs to local or window-based approaches.

B. General Scheme

GeFolki takes two input images resampled at the same rate, and outputs the two components of the flow at every pixel position. The aim of the algorithm is to minimize a given criterion J(u; x) [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF]:

J(u; x) = x ∈S ω(x -x)(f 1 (I 1 (x )) -f 2 (I 2 (x + u(x))) 2 (3)
where ω defines a local window of radius r and size (2r + 1)x(2r + 1), thus ω(x) = 1 iff |x| ∞ ≤ r and 0 else; f 1 is a function applied to the master image I 1 , and f 2 is a function applied to the slave image I 2 . These functions f 1 and f 2 are designed to project the images in a space where f 1 (I 1 ) and f 2 (I 2 ) are similar enough to validate the brightness constancy model. They significantly improve the robustness, but are also designed to retain the efficiency of our method.. Multi-scale algorithms can deal with displacements larger than a few pixels because they initialize the flow by estimates from coarser image scales, where displacements are small enough to be estimated by local optimization. Unfortunately, too much down-sampling can remove information that may be vital for establishing the correct matches. Therefore, to converge to a correct solution, the initialization is important. We propose to start with a couple of images first computed in a common resolution grid. We compute the projective transformation based on the selection of four tie-points. The outcome from this first coarse registration becomes our input images for registration.

The first step is to create a pyramidal structure of downsampled images, in order to be able to find large displacements. The number of pyramid levels L affects the amplitude of the desired movements. Then, we have to solve the minimization of the criterion J(u; x) iteratively in each scale, as represented in Fig. 3. The solutions are incrementally refined over the scales. Then, for each pyramid level, from the coarsest resolution to the finer:

• first, functions f 1 and f 2 are applied to the intensity images I 1 and I 2 . This functions are detailed below; • then, the algorithm minimizes the criterion expressed in eq. 3 by first-order Taylor expansion of the residuals at the current flow estimate. At the first pyramid level, this initial current flow value is 0, and at the other levels, it derives from the expansion of a previous low resolution estimate. The quadratic term of eq.3 is a composition of several functions with non-trivial derivatives, for example the rank function that is not continuous. Recall, that the Gauss-Newton iterations require a linearization around the current solution.

Here, we cannot compute these derivatives explicitly. Instead we approximate the derivatives numerically, with regard to the current discretization of the image function. In other words, we directly consider f 1 (I 1 ) and f 2 (I 2 ) as functions defined on the image grid, after applying the similarity correcting transformations f 1 and f 2 , and operate with linear finite elements defined on the image grid. We did not experience any problems with this simplification.

In GeFolki, the function f 1 and f 2 are designed to ensure the robustness of the algorithm, as well as the execution speed. The success of the criterion lies in the presence of macroscopic components similar in both images. It is imperative to adjust the parameters in such a way that these structures are not affected by noise or by smaller details that differ in texture. Thus, in GeFolki, f 1 and f 2 are defined as follows:

• For coregistration of SAR images, function

f 1 = f 2 = R
is a rank filter, as proposed in [START_REF] Plyer | A New Coregistration Algorithm for Recent Applications on Urban SAR Images[END_REF]. This function will be explained below. Note that in this case, we are brought back to the initial eFolki algorithm, evaluated in [START_REF] Plyer | A New Coregistration Algorithm for Recent Applications on Urban SAR Images[END_REF]. • For coregistration of heterogeneous images, f 1 and f 2 will include, in addition to the rank filter:

a filtering step based on the rolling guidance filter, in order to make image textures more similar, a function of local contrast inversion, for all cases where the contrast exhibited between the different sensors is not the same. It is very often the case for co-registration between optical and radar images, and more generally for images acquired in different parts of the electromagnetic spectrum, for which the scattering levels change, as shown in Fig. 5. Therefore, in this case f 1 and f 2 are not strictly the same, and can be written as:

f 1 = R • g f 2 = C • R • g
where g is the rolling guidance filter, R is the rank-filter, and C is the local contrast inversion decision function.

We choose for all of them fast execution functions to maintain the speed. They are described in the following sections.

C. The Rolling guidance filter g

An issue concerns the difference in texture and resolution in the original images. When the resolution of an image is too high compared to the other one, then the textures are generally different, and this makes the coregistration more difficult to do. Moreover, a detail present in one image, and not visible in the other one, can perturb the result of flow calculation, as well as the presence of speckle in SAR images. For example in the case of registration of SAR-LIDAR images, the main difficulty lies in the differences existing about texture and resolution effects. These differences make the flow calculation setting particularly challenging.

For these reasons, in order to maximize the efficiency of the algorithm, we propose to apply a pre-filtering step that can make the intensity images more comparable. Among them, we have successfully tested the Rolling Guidance Filter [START_REF] Zhang | Rolling guidance filter[END_REF], that removes small-scale structures but maintains the edges in the image. To that end the image is first blurred with a Gaussian filter. This step essentially removes small structures and artefacts. Starting from this blurred version, repeated bilateral filtering with the original image itself, then recovers the edge information.

In our experience, Rolling Guidance Filter enables to filter images with the complete control of detail smoothing under a scale measure, which is especially interesting to use in GeFolki since it is also a multi-scale strategy. Finally, the Rolling guidance filter implementation achieves real-time performance and then can be integrated to our algorithm without degrading its computing performance. The complete description of the algorithm is given by [START_REF] Zhang | Rolling guidance filter[END_REF], and the source code is provided by the authors.

D. The Rank filter R

This rank function R(I) applied to the image I, is expressed as:

R(I)(x) = #{x : x ∈ S R (x) with |I(x)| > |I(x )|}, (4)
where S R (x) is a neighboorhood of the pixel x.

We use # to denote the number of elements of the set. This rank function consists in replacing the intensity value I(x) by the number of pixels in the neighborhood of x with a intensity lower than I(x). Thus, the rank transform depends on the relative order of pixel values within a local neighborhood, which is invariant under various transformations of the data, namely those which preserve this order. The effect of the rank transform is a non-linear filter that highly compresses the signal dynamics. By this compression effect on the signals gradient, the rank filter enhances the robustness of the motion estimation. This has been already studied in [START_REF] Zabih | Non-parametric local transforms for computing visual correspondence[END_REF], and used in [START_REF] Eklundh | Dense optical flow by iterative local window registration[END_REF]. Other functions such as Census transform could also be used, but empirically deliver worse results. Moreover, rank filter is better suited to a massively parallel implementation.

E. The Local Contrast Inversion decision function C

In order to be compliant with the brightness constancy model, we need to address the possible differences in terms of contrasts. From one image to another, these contrasts can be locally reversed. This is particularly the case between a radar image and an optical image.

In order to decide whether the intensity value of the pixel under consideration x 0 should be inverted or not, we proceed as follows. We consider W a (2M+1)×(2M+1) sliding window centered around x 0 , defined by

W(x 0 ) = {x : |x -x 0 | ∞ ≤ M}.
We determine whether the set of points {(I 1 (x), I 2 (x + u(x))) : x ∈ W(x 0 )} is closer to the line defined by the equation y = x or the line defined by the equation y = 1 -x. To this aim, we have to compute the following criteria:

C 1 (x 0 ) = x i ∈W(x0) |I 1 (x i ) -I 2 (x i + u(x i ))| (5) C 2 (x 0 ) = x i ∈W(x0) |1 -I 1 (x i ) -I 2 (x i + u(x i ))| (6)
If C 1 (x 0 ) < C 2 (x 0 ) we do not have to invert the intensity value of the second image, whereas if

C 2 (x 0 ) < C 1 (x 0 ) we replace I 2 (x 0 + u(x 0 )) by 1 -I 2 (x 0 + u(x 0 )) in image I 2 .
One example for each case is given on Fig. 4 and 5 with M=100. The first case concerns images of a transition between a vegetated area and a water one. In both radar and optical images, the contrasts are the same, and the intensity point cloud is closer to the line y = x than the line y = 1 -x. As a consequence, no inversion will be applied to the central pixel. In the second case, we see a transition between a forest area and a bare soil area. In the optical image, forest is darker than the soil, whereas forest is brighter than the soil in the radar image. Thus, the corresponding point cloud is closer to the line y = 1 -x. The value of the central pixel x 0 of the slave image will be replaced by 1 -I 2 (x 0 + u(x 0 )).

One situation where this strategy does not work could be when one of the images contains an element, and the other image does not contain it. The strategy of parameter settings for the window size will be discussed in section IV-A. For the whole images, these criteria can be expressed as:

C 1 (x) = (|I 1 -I 2 | * W)(x), (7)
C 2 (x) = (|1 -I 1 -I 2 | * W)(x) (8)
where * denotes the convolution between images and W is a (2M+1)×(2M+1) all-ones matrix. The advantage of these criteria is that they can be rapidly computed by using integral images.

Then the transformed slave image C(I 2 ) is computed from:

C(I 2 ) = (C2 < C1) • I 2 + (C1 > C2) • (1 -I 2 ) (9)
where (C2 < C1) and (C1 > C2) are simple logical arrays corresponding to the respective conditions expressed for each pixel, and • is the Hadamard product.

III. RESULTS

A. Description of the data set

Now that we have presented the algorithm, we apply it to a concrete case of interest: the heterogeneous image registration on the site of Remningstorp, a Swedish forest. Firstly, this site was chosen because the forest is a particularly interesting case for fusion; and many open dataset exist for the scientific community. Secondly, it represents an interesting case for coregistration because forest images contain few specific features such as segments and isolated points. Finally, the images are airborne high resolution images, where deformations are particularly interesting to estimate with a dense method. The data set contains a radar image, an airborne Laser Scanning (ALS) image, as well as an optical image.

The SAR (Synthetic Aperture Radar) acquisition is a P-band full polarimetric image taken by the airborne Onera system SETHI [START_REF] Bonin | The Airborne SAR-System: SETHI -Airborne Microwave Remote Sensing Imaging System[END_REF] on 23th September 2010 during the BioSAR 2010 campaign. It corresponds to the first track of the BioSAR delivery. This acquisition has been chosen to maximize the common footprint with the LIDAR data provided by SLU (Swedish University of Agricultural Sciences). The quality of localization has been estimated on a trihedral corner with an accuracy of about one meter. Accuracy of trajectory is 20cm. Range and azimuth resolutions are 80 cm. Local incidence angle varies from 25 to 60 degrees, and the image is 12,250 x 7,000 pixels.

The LIDAR acquisition was carried out by the Helicoptermounted TopEye Mk III system on 29th August 2010. All LIDAR data have been processed by SLU, i.e. going from the collected LIDAR raw data of point clouds to various products delivered as several raster files. Here we use a raster image with a 0.5mx0.5 m resolution, which contains the difference of height value of the first LIDAR return and height value of the DEM (Digital Elevation Model) at the LIDAR return horizontal position. This image is called DSM (for Digital Surface Model) image and is georeferenced in WGS84.

For the georeferencing of the SAR data, another DEM was provided by FOI, the Swedish Defence Research Agency. This DEM covers the entire radar images. At the date of the BIOSAR 2010 campaign, this DEM was a commercial product with a single user license and not included in the deliverable distributed by ESA for public access. However, since July 2015 it is provided without license with simplified terms of use by the National Land Survey as open data. It contains ground height values at the intersections of a 50 m grid. This product is referred as GSD50+ and is described in [START_REF]GSD-Elevation data, grid 50+ hdb[END_REF]. Standard error in height is announced to be 2.5 m. To test the influence of this DEM on the coregistration, in particular on georeferencing, we have also used the product delivered by SRTM90 data before 2015, which gives ground height values on a 90 m grid.

B. Tested co-registration scenarios

For the coregistration of two images, we have to choose one as master and the other one as slave. When a SAR image is included in the data set, we choose this one as the master image for the whole coregistration products for at least two reasons:

• Projecting a radar image on the ground introduces artefacts, because of the special SAR geometric configuration [START_REF] Pétillot | Radar-coding and geocoding lookup tables for the fusion of GIS and SAR data in mountain areas[END_REF]. These artifacts are further amplified for high relief terrains, and also when the incidence angle variation is great, such as encountered for our airborne acquisitions. • The choice of a master SAR image is in favor of advanced processing methods such as interferometry that will not be affected by a resampling of the radar image.

In order to enhance the usefulness of GeFolki for the coregistration, we have tested several scenarios as follows:

• The SETHI SLC data have been delivered by ONERA with an auxiliary transformation grid file to be able to link precisely one pixel on the image to one geographic location. This grid file allows to link a pixel position with a height to a precise latitude and longitude. It can also be used to associate an image position and a height to a geographic position. Therefore, we have applied the process to link each pixel of our SAR image to a geographic location, by using different available DEM. This process is described in [START_REF] Ulander | BIOSAR 2010 Technical Assistance for the Development of Airborne SAR and Geophysical Measurements during the BioSAR[END_REF]. Then, this geocoding information is used in combination with the LIDAR one to find the differences of positions of image pixels in map coordinates. Thereafter, it is possible to resample one image onto the other one. Results obtained by this process will be called coregistered products by geocoding.

The corresponding synoptic is illustrated on Fig. 6. This method does not utilize the flow algorithm. It has been implemented in order to compare the results with those obtained using GeFolki. It also allows to judge the influence of the quality of the DEM on the accuracy of registration. Indeed, we have applied it using the two types of DEM available on the whole image: SRTM90 and GSD50+ provided by the National Land Survey. • In parallel, we register images without using any georeferencing data. We begin with images after the initialization procedure explained in section II B., as represented in Fig. 7. Then we apply GeFolki to compute the flow, and the slave image is re-sampled accordingly. We call this approach GeFolki after coregistration by manual initialization. This approach has been chosen because it enables us to attempt a coregistration even if no auxiliary data is available.

As a baseline we use a method whose similarity metric is based on mutual information [START_REF] Kroon | Mri modalitiy transformation in demon registration[END_REF]. Because the baseline algorithm is more time consuming, we could not apply it at full resolution (12,250 x 7,000 pixel). For that reason, we have restricted the test on a part of the image only, and thereby, only to qualitative evaluation. The initial algorithm introduced by Thirion [START_REF] Thirion | Image matching as a diffusion process: an analogy with maxwell's demons[END_REF], called demon registration, is a popular algorithm for non-rigid image registration because of its linear computational complexity and ease of implementation. It is a non parametric registration that approximately solves the registration problem by successively estimating force vectors that drive the deformation toward alignment, and smoothing the force vectors by Gaussian convolution. In this version, the method has been adapted to heterogeneous images. This algorithm is provided for free by the authors in matlab language. Therefore, it is close in spirit to GeFolki, in the sense that it handles the conventional formulations of optical flow problem, adapts the similarity criterion to heterogeneous images, and is able to compute the flow for each pixel in the images.

The conclusions drawn from our various tests are:

• Geocoding methods with our data set are insufficient to ensure a pixel precision. The impact of DEM accuracy on the result is clearly demonstrated in next subsection. • The improvement brought by GeFolki is significant. It can be particularly useful for all applications where the highest accuracy is required, for instance for the counting of trees. • GeFolki was sucessful even with its coarse initialization trough an projective transformation based on a few tiepoints.

• With the same coarse initialization, the baseline algorithm based on mutual information was not successful. Therefore, we have managed to use it in combination with our best geocoding result as initialization. The results appear to be still of inferior quality and needs a 20 times more computation time. In next subsections, we detail these qualitative results, first, of the LIDAR DSM image at the best resolution (0.5m) on the SAR image, and then of an optical aerial photography on the SAR image.

C. 0.5mx0.5m elevation map with the SAR image

We consider here coregistration results of a 1,000 x 1,000 pixel part of the whole image SAR image, with the LIDAR 0.5mx0.5m elevation map. This part contains forest stands, isolated trees, electrical cables and pylons. Each result is presented as a colored composition in which the radar image is in the green channel and the coregistered LIDAR image in the magenta channel. When a pixel is white or grey, we have similar intensities in both images. If not, a pixel will tend to take the color of the channel with the highest intensity value. With this representation, it is then possible to appreciate the quality of a registration, in particular when we focus on isolated scatterers.

The two first results presented on Fig. 8 are coregistered products by geocoding. For the first one, we have used the Figure 6. First coregistration scheme by geocoding from a geometric point of view: the latitude and longitudes are computed for each pixel of the images according to topographic and track information. Then the flow is computed between these both coordinates, and the slave image is resampled accordingly. SRTM90 DEM, and for the second one the GSD50+ DEM which is more accurate. For both cases, we can notice that a residual shift remains, notably between forest edges and small structures such as pylons. These results also demonstrate that the precision obtained by geocoding is better with the most accurate DEM of the Land National Survey, than with the SRTM product. Thus, the DEM choice is crucial on the final accuracy.

(1)

(2) Figure 8.

Result of coregistation for DSM LIDAR data on a 1,000 x 1,000 pixel SAR image sample over Remningstorp forest (Sweden) (1) after coregistration by geocoding using SRTM, (2) after coregistration by geocoding using GSD50+ DEM Then, Fig. 9 presents the results obtained by the baseline algorithm based on mutual information and by GeFolki. The baseline algorithm has been initialized by the previous result of best geocoding, because it was diverging for the coarse initialization based on the projective transformation.

For these two new results, the images show less saturated colors, which means that the coregistration results are better than using the geocoding method. However, some colored pixels remains, mainly because of the lateral viewing geometry of the radar, that implies that the return of the tree is often clearly double: the branches echoes appear in HV polarization at the first range cells, whereas the trunks double bounce echo appears in HH-VV polarization at further range cells. Due to its vertical illumination, the LIDAR return of one tree is a single patch. These effects are illustrated on Fig. 10. Note that for this reason, we perform the computation of the flow by using the radar image calculated in HH-VV polarization, in order to make the double bounce echo of the ground matching to the LIDAR return echo.

(1) (2)
Concerning the baseline algorithm based on mutual infor- 11 show that tree superposition is good for both GeFolki and baseline algorithm; however cables are not straight any more in the case of the coregistration by baseline algorithm.

The execution time of the baseline method in matlab is 248 seconds, while Gefolki needed only 7 seconds for the same images (1,000 x 1,000 pixels).

D. Optical and SAR images

Our algorithm GeFolki is now applied for the registration of radar and optical images. The optical image used in this test is distributed by Lantmäteriet. Several areas of the scene present intensity contrasts that differ between the SAR image and the optical one, as we will show later in Fig. 16. Forests are dark patches in the optical one, and bright in the radar images. Water surfaces are dark in both cases. The intensity of the soil varies greatly in both images but not necessarily in the same way for the radar and for the optical sensors. In this case, success of GeFolki lies mainly in the inclusion of the local contrast inversion.

Optical and SAR images now implies even more different geometries. The optical scene is viewed from a given camera position where the main source of light is the sun. In the case of our images, the camera is almost at the top of the tree, the sun is low in the sky and consequently the shadows are relatively large. The optical signature of a tree will consists of the crown, which partially hides the foot of the trunk, and the shadow. The crown has little contrast with the ground; it is especially the shadow that is visible in the case of isolated trees, as seen as in Fig. 12 on the top.

Also, in order to create a visualization of the differences between both radar and optical images, we first reverse the contrast of the optical image before assigning the magenta color channel in the output image. The green color is assigned to the radar image in the final composite visualization. In the case of a perfect coregistration, the colored overlay makes us see the optical shadow in magenta, whose direction crosses the radar double bounce echo, colored green. The choice of the colors is explained in Fig. 12.

Using the same procedure as in the case of LIDAR to coregister the optical image with the SAR image, we obtain the different results presented in Fig. 13 for geocoding, baseline registration and GeFolki.

Once again, the coregistration algorithm applied to remote sensing heterogeneous images results in an interesting product in terms of spatial accuracy. The algorithm based on mutual information also improves the coregistration results in most areas, albeit unevenly, as shown in Fig. 14. On the first extract, both baseline registration and GeFolki improve the matching of trees, whereas in the second case, only GeFolki

IV. ROBUSTNESS AND PARAMETER SETTING ANALYSIS

Several parameters have to be set for the whole coregistration process: the window size parameter for the local inversion contrast, and the parameters used for the calculation of the flow. We study these two kinds of parameters independently, to evaluate their impact on the quality of results.

A. Discussion about the local contrast inversion

The first point is that for the coregistration of optical and radar images, the algorithm fails to estimate the flow without the local contrast inversion. Indeed, estimation of the flow without the contrast inversion function diverges, even if a good initialization is done. Therefore, this function is crucial in the algorithm.

(1)

(2)

(3) The criterion for inversion depends on the size of a local window. Therefore, we must keep in mind that this size is the outcome of a compromise between:

• the scale of structures offering contrasts that we want to correct: the size of the window has to be in the same order than the size of these structures; • the flow to estimate: when a significant displacement between the two images exists, the inversion process must be performed at a larger scale than this displacement, otherwise the algorithm may fail on the transition areas. Ultimately, through several tests on different images, we can confirm what we expect intuitively: similar to the flow calculation, the contrast inversion is very effective when conducted in a multi-scale manner. At a given pyramid level, and for window size equal to that of the local window ω in eq. 3, the contrast inversion is conducted at a scale large enough compared to the flow, and therefore is always carried out with an appropriate setting.

To illustrate this, we take a sample of our optical and radar images. The initial image pair is presented in Fig. 16, simply in native gray scale formats.

At the bottom, we show intermediate results of contrast inversion, obtained for images that have been scaled down by a scale factor 8 through a dyadic pyramid reduction of three levels (2 3 = 8). The images have been successfully transformed by the local contrast inversion, for a radius equal to 8, that corresponds to a 17x17 local window (2r+1 x 2r+1). Globally, the obtained images are very similar in terms of intensity levels. Some small areas are still different in contrast, such as the small lake highlighted by red rectangles. Nevertheless, these images depict the solution at an intermediate stage and when we apply the inversion algorithm to images at the initial resolution, we will see that the result is effective on these areas, as shown on Fig. 17.

In conclusion, a multi-scale strategy is essential for the contrast inversion decision, because this one is robust only for local window large enough, and displacements small enough A -Image display with mosaic after coregistration by geocoding B -Image display with mosaic after coregistration by GeFolki

B. GeFolki parameter settings

Now that the sensibility of the contrast inversion step has been analyzed, we want to discuss the influence of the parameters of GeFolki in the case of remote sensing heterogeneous image coregistration. The main features are:

• r, called the radius, that describes the size of the local window; • L, the number of pyramid levels;

• K, the number of iterations. Empirically, the found solutions were not very sensitive to K. A prior sensitivity analysis was conducted in the case of SAR images in [START_REF] Plyer | A New Coregistration Algorithm for Recent Applications on Urban SAR Images[END_REF]. It was shown in particular that it was preferable to use large radii in order to minimize the impact of the speckle, this latter being a level of texture that provides no useful information for the correlation between images.

However, the SAR images were acquired in very similar conditions. The flow does not vary much locally. In the case where the images are taken with different geometries, the terrain projection effects can induce more local deformation between the two images, and in this case, it could be necessary to use a smaller radius to be more precise.

Still in the case of SAR images in close acquisition configurations, the desired displacements are small. In this case, it is sufficient to reduce the number of pyramid levels to only two. In the other cases investigated in this paper such as heterogeneous images, we encountered situations where even the coarse coregistration in first stage leads to remaining errors of several tens of pixels to be compensated after. In such cases, more pyramid levels are needed to estimate these major shifts.

However, these two parameters, radius and number of levels, are intrinsically linked. If the overall displacement to find is large, a high number of pyramid levels is required. Moreover, a larger window leads to improved matching robustness but also constrains the motion which can be estimated, what would not be appropriate especially at a coarse resolution. This confirms that the quality of initialization plays a significant role in the successful outcome of GeFolki.

Overall, we therefore recommend the following procedure in order to optimize the effectiveness of our flow algorithm:

• The initial pyramid level is chosen directly based on the expected residual displacement after initialization. For example, if the flow to estimate remains greater than several tens of pixels, we have to choose a rather high pyramid level number, typically equal to 6, equivalent to a 1:64 scale. On the contrary, if images are already very close after the initialization step, the pyramid level can be only 2.

• to choose relatively a rather large matching window. An effective strategy is to explore several values to ensure robustness, for example a list of 6 radii ranking from 4 to 24 by step of 4.

V. COREGISTRATION PERFORMANCE

To date, there is no universal method to quantify the accuracy of a registration algorithm for heterogeneous images. Also, we propose a methodology for quantitative assessment of performance. The key concept is to measure distances between paired features in the two registered images. These features will be different according to the image pairs under consideration.

A. Metrics for LIDAR-SAR coregistration

The main idea is to detect trees in both images, to match them, and to to measure the distance in pixels for each defined pair. This pixel-wise distance can be converted into physical distance. Depending on the number of points considered, it is also possible to carry out a statistical analysis of the results. Given the nature of our images, we propose to calculate performance by considering the positioning on the tree bases on the ground.

In practice, we proceed in three steps, described in Fig. 18:

Figure 18.

The different steps to evaluate metrics on LIDAR-RADAR coregistration results

• First in the LIDAR image, we extract the individual positioning of each tree. This positioning of individual trees from the LIDAR DSM is performed using the package rLiDAR [START_REF] Silva | Package rlidar[END_REF]. The two processing steps are a mean smoothing in order to eliminate spurious local maxima caused by tree branches then non-maximum suppression within two meters distance. Both functions have been used with a fixed 5x5 windows. • We extract in the polarimetric radar image the isolated double bounce echoes. This detection is carried out starting from intensity image polarization HH-VV. A first simple threshold detection is performed on this image. The binary mask thus created is dilated by using a 3 radius disk structuring element object. We then consider all connected components found after this dilatation. We keep among these components, only those whose number of pixels does not exceed 100. Otherwise, we suspect several trees to be aggregated. For each detection area, the coordinates are calculated at the maximum intensity pixel.

If several adjacent pixels have the same peak intensity, we take the center of gravity of these pixels.

• The last step is to search for each trunk detected in the master image, position of its correspondent, supposed the nearest in the second image. Finally we calculate the distance between the master pixel position and its corresponding position. The LIDAR detection performance seems very good, while radar detection fails once the trunk density is high. For this reason, we have manually selected areas where detection results are good. It is also possible to automatically select the areas where the tree densities found in the LIDAR data are low. Five areas selected are indicated on Fig. 19 on the radar image on the right. From figures 20 to 24, we present the results relative to these areas. For each of them, we compare the results in terms of distance, after a coregistration performed:

• by geocoding following the scheme described in Fig. 6.

• after a simple projective transformation based on four Ground Control Points, as the first step of the scheme described in Fig. 7. This result constitutes the initialization of GeFolki. • after coregistration of this latter product by GeFolki, as the final product described in Fig. 7. Each image is as follows:

• The first line shows the colored compositions of the LIDAR and radar images: one codes the magenta channel, and the other the green channel. This representation allows to visualize the effect of any mismatch between images. • The second line for figures 20 and 21 indicates the trees location of the LIDAR image (red circles) and the radar image (black dots). • The last line is a table where we summarize the mean distance, the standard deviation, and the root mean square error obtained for all detected trunks. Thus, Fig. 20 concerns a simple case of five well isolated trees whose LIDAR or radar signatures correspond to easily identifiable echoes. Clearly on the colored composition, we see that GeFolki allows to match the radar echoes with the crown that can be seen in the LIDAR image. Quantitatively, the geocoding method brings up a residual gap between the two images of about 10 pixels. On the image obtained by projective transformation, the calculated distances are of the same order of magnitude. However, in this case it is very likely that the error is largely underestimated: the found matches are not correct, because the gaps between trees are less than the error of existing offset in this area. Finally, GeFolki leads to a good match with a root mean square error around one pixel. The error is below the expected accuracy of our trunk detection chain. On the second example represented on Fig. 21, we see a row of trees arranged in pairs. The order of magnitude of error is the same than previously.

However, we can see in this case that our evaluation can be more pessimistic for a precise algorithm than the real quality of coregistration: in fact, some trees are not detected in the radar image. For these missing trees, the distance measured on the LIDAR image will be done with the nearest tree, but not the true corresponding one. This will tend to increase the average distance calculated for GeFolki. On the contrary, this matching error has the opposite effect for the projective case: in this case, the nearest tree is not the real one and the metrics leads to underestimate the error, as in Fig. 21.

Fig. 22 is an example chosen along power lines. On the left there is one row of trees, and on the right side we see a group of trees. The evaluation leads to a slightly worse score for GeFolki than in the previous examples, but still better than the other two methods. Visually, it is mostly trees arranged on the left that make them slightly less good statistics.

Another example of grove is given in Fig. 23. It is located in a part of the image very far from the previous one, at a distance of about 5 km. In this area, the projective transformation fails completely. Geocoding and GeFolki have equivalent statistics, slightly in favor of the latter, in terms of root mean square error about 8 pixels for Geocoding and 6 pixels for GeFolki. Note that in this example, the displacement obtained after the simple projective model is rather strong, and GeFolki can still recover from such an inaccurate initialization.

A final example with more trees is shown in Fig. 24, still in favor of the GeFolki method.

B. Metrics for Optics-SAR coregistration

For the optical images, the detection of trees is much more difficult to handle. For this reason, we first go through a manual selection of some reference points, at the foot of the various pylons detectable in the images. Most of them are found in pairs. In this case, the coordinates of each of them are then averaged. We then get 7 coordinates, highlighted by yellow rectangles in Fig. 25, each corresponding to a pair of pylons. Note that in the case of radar, the echo at the foot of the tower is very strong, as can be seen in Fig. 26. Furthermore, the selection of the point is further facilitated by the polarimetric color signing the double bounce mechanism. Contrariwise, the position of the pylon base is more uncertain in the optics image, because these objects provide little contrast with the ground. Despite these difficulties, we are able to locate these targets which are well distributed throughout the image extent; furthermore echoes are sufficiently focused.

Then, we calculate for each pylon the difference between the coordinates of the radar and the optical image. This procedure is applied on the optical image coregistered by geocoding, and the one coregistered by GeFolki. The results are given in Table 27. We can note that GeFolki presents the best results.

Finally, we also looked at the errors of coregistration on other types of objects: the outline of the lakes and some rows of trees. We performed manual reading distances in the 6 areas circled in yellow on the mosaic of Fig. 15. Four distances are recorded at the transition on rows of trees, and two distances are measured on the edges of the lake. These distances lead to a root mean square error of 18.5 pixels in the case of data coregistered by geocoding, and 0.8 pixels in the case of data coregistered by GeFolki.

In summary, the GeFolki method provides more accurate results than geocoding, with pixel-precision, compared to the ten pixels from our method of geocoding.

VI. CONCLUSION

This article introduces the use of an optical flow algorithm called GeFolki to register heterogeneous remote sensing images. This consists in a dense flow method derived from the Lucas-Kanade algorithm, with a multi-scale strategy, applied to images transformed by a rank-filtering, and, for heterogeneous images, also by a rolling guidance filter and a contrast inversion decision method.

The algorithm has been applied and qualitatively compared to coregistration made using only geocoding information, and also to an alternative method of baseline registration based on • more accurate than geocoding on all our examples;

• free from some defects such as deformation induced by the baseline registration;

• more robust than this same method with respect to the initialization. Quantitative performances have been derived, both for a coregistration based on geocoding, and for GeFolki. Performances are based on detection of trees for the SAR/LIDAR coregistration, and on manual measurements on some pylons, tree rows and lake edges for the SAR/optical coregistration.

The registration precision of GeFolki is about a few pixels, (1 pixel=80 cm), compared to ten of pixels by geocoding. This precision opens new opportunities for studies related to data fusion, for example, with applications involving the counting of trees. It appears that the fusion of different images may improve detection performance, but only if we can trust the accuracy of registration.

Finally, GeFolki has been proven to be faster than the baseline registration method, by a factor of 10 to 20.

Our next efforts will focus on:

• Taking into account the several bands of multivariate images, as opposed to a panchromatic image or single polarization radar image. Indeed, Hyperspectral or LIDAR data contains several bands having different contrasts, as well as the polarimetric radar images with the different polarimetric channels. To date, the choice of the band has not been reasoned. But a work could be undertaken on how to make the most of this data diversity. • The possible extension to more complex scenarios, for example with two images containing mountainous terrains or buildings, and acquired with different headings. This case is particularly difficult because these objects will be projected on the ground differently following the acquisition mode.

VII. ACKNOWLEDGMENT All data used in this paper come from various institutions:

Figure 1 .

 1 Figure 1.Definition of variables used in our optical-flow algorithm: To register the slave image I 2 to the master image I 1 , we model the displacements between both images u(x) for each pixel x in the master image. Corresponding pixels I 1 (x) and I 2 (x + u(x)) are supposed to be close with respect to a similarity measure.

Figure 3 .

 3 Figure 3. To enable the recovery of large displacements, a multi scale strategy is employed: After the generation of image pyramids for master and slave image, the optical flow is computed from coarse to fine, by minimizing the cost function (eq 3).

Figure 4 .

 4 Figure 4. Example of the contrast inversion decision. Left: Optical master image I 1 for a local window. Middle: Radar slave image I 2 for the same window. Right: Scatter Plot of I 1 (x i ), I 2 (x i ) for the example. Our criterion defined in eqs. 5 and 6 leaves the value of the central pixel x 0 of the window unchanged.

Figure 5 .

 5 Figure 5. Example of the contrast inversion decision. Left: Optical master image I 1 for a local window. Middle: Radar slave image I 2 for the same window. Right: Scatter Plot of I 1 (x i ), I 2 (x i ) for the example. Here, our criterion defined in eqs. 5 and 6 inverts the value of the central pixel x 0 of the investigated window.

Figure 7 .

 7 Figure 7. Coregistration scheme by optical flow calculation: first a coarse coregistration is performed by warping the slave image on the master image trough a projective transformation based on the manual selection of four tie-points. This first step constitutes the initialization of GeFolki. Then, GeFolki computes the flow for this couple of images, and the final slave image is resampled accordingly.

Figure 9 .

 9 Figure 9. Result of coregistation for DSM LIDAR data on a 1,000 x 1,000 pixel SAR image over Remningstorp forest (Sweden) (1) after coregistration by an algorithm based on Mutual Information, (2) after coregistration by GeFolki

Figure 10 .

 10 Figure 10. Effects of geometry on the signature of a tree in the scene for polarimetric radar images and for LIDAR image. Example on a very small 30 x 40 pixel extract of the different products containing four trees: polarimetric SAR image on the left, DSM LIDAR image on the middle, and superposition after our best coregistred product on the right

Figure 11 .

 11 Figure 11. Result of coregistation for LIDAR data and SAR image over Remningstorp forest (Sweden) for two 100 x 100 pixel extracts (1) after coregistration by geocoding (2) after coregistration by the baseline algorithm based on Mutual Information, (3) after coregistration by GeFolki

Figure 12 .

 12 Figure 12. Effects of geometry on the signature of a tree in the scene for the optical image. Example on a very small 50x50 pixel extract of the different products containing some trees: initial optical image on the left, then inverted optical image, SAR image, and superposition after our best coregistred product on the right

Figure 13 .

 13 Figure 13. Result of coregistation of an optical image on a 1,000 x 1,000 pixel SAR image over Remningstorp forest (Sweden) (1) after coregistration by geocoding (2) after baseline coregistration based on Mutual Information, (3) after coregistration by GeFolki

Figure 14 .

 14 Figure 14. Result of coregistation for optical data and SAR image over Remningstorp forest (Sweden) for two 100 x 100 pixel samples (1) after coregistration by geocoding (2) after coregistration by an algorithm based on Mutual Information, (3) after coregistration by GeFolki

Figure 15 . 1 ,

 151 Figure 15. 1,600 x 1,500 pixel mosaic of the polarimetric SAR image and the RGB aerial photography. A: by geocoding; B: by GeFolki

(Figure 16 .

 16 Figure 16. Top left (A): the Red Component of the optical image -Top right (B): the SAR image in HH channel. Bottom left (C): the Red Component of the optical image after local inversion procedure -Bottom right (D): the SAR image in HH channel after contrast enhancement.

Figure 17 .

 17 Figure 17. Zoom of contrast inversion on images at initial resolution

Figure 19 .

 19 Figure 19. Selection of areas for metrics in coregistration of LIDAR and SAR images-Left: the LIDAR image, Right: the 12,250 x 7000 pixel polarimetric image with areas selected for our metric highlighted by yellow rectangles

Figure 20 .

 20 Figure 20. Performance evaluation for area 1

Figure 24 .

 24 Figure 24. Performance evaluation for area 5

Figure 25 .Figure 26 .

 2526 Figure 25. Selection of areas for metrics in coregistration of optical and radar images-Left: the polarimetric image, Right: the optical image with seven pairs of pylons highlighted by yellow rectangles

Figure 27 .

 27 Figure 27. Examples of flow errors at reference points

  • LIDAR data have been provided by Johan Fransson from SLU, Swedish University of Agricultural Sciences, and Lars Ulander, from FOI, Swedish Defence Research Agency, for whom we are deeply grateful. All LIDAR images have been processed by SLU. • Polarimetric SAR data set over Remningstorp is provided from the BIOSAR 2010 campaign funding by ESA, acquired by Onera. • The optical image over Sweden has been provided by Lantmäteriet. One of the authors, Guillaume Brigot, is a PhD student jointly supervised by Elise Koeniguer in Onera and Marc Simard in JPL. His PhD is supported jointly by TOTAL and Onera. The authors are thankful to Malcolm Davidson and Michael Foumelis from ESA-ESRIN for their very informative recommendations. We are also immensely grateful to Guy Le Besnerais for his comments on an earlier version of the manuscript, although any errors are our own. Finally, we warmly thank Charlotte Smetanka for her careful reading of the manuscript.


Yu-Liang Liu
email: au4377@au.edu.tw

Routing and Wavelength Assignment for Exchanged Hypercubes in Linear Array Optical Networks

Keywords: WDM optical networks, Routing and wavelength assignment, Exchanged hypercube, Linear array, Congestion

The exchanged hypercube, denoted by EH(s, t), is a new interconnection network obtained by systematically removing links from the hypercube, while preserves many appealing properties. This paper addresses the routing and wavelength assignment for realizing exchanged hypercubes communication patterns on linear array WDM optical networks. By using congestion estimation, we derive a lower bound of the minimum number of required wavelengths, and propose an optimal wavelength assignment algorithm that uses 2 s+t-1 + 2 t /3 wavelengths.

Introduction

In a wavelength division multiplexing (WDM for short) optical network, the bandwidth in optical fiber is partitioned into multiple virtual channels, in which different stream of data can be transmitted simultaneously using separate virtual channels. In this context, a virtual channel corresponds to a wavelength. In general, a WDM optical network consists of routing nodes interconnected by point-to-point fiber links. To achieve all-optical communication without optoelectrical conversions at intermediate nodes, end-to-end lightpaths are usually set up between each pair of source-destination nodes. A connection or a lightpath in a WDM optical network is an ordered pair of nodes (S, D) corresponding to transmission of a packet from source node S to destination node D.

The primary issue for WDM optical networks is to select a proper path and wavelength satisfying the wavelength-continuity constraint and the distinct wavelength constraint for each connection of a given communication pattern so that the number of used wavelengths is minimized [START_REF] Yu | Routing and wavelength assignment for 3-aryn-cube in array-based optical network[END_REF][START_REF] Yu | Routing and wavelength assignment for 3-ary n-cube communication patterns in linear array optical networks for n communication rounds[END_REF][START_REF] Yuan | Optimal routing and channel assignments for hypercube communication on optical mesh-like processor arrays[END_REF][START_REF] Zang | A review of routing and wavelength assignment approaches for wavelength-routed optical networks[END_REF]. Up to now, there have been some works about routing and wavelength assignments in optical networks [START_REF] Chen | Routing and wavelength assignment for hypercube in array-based WDM optical networks[END_REF][START_REF] Chen | Wavelength assignment for realizing parallel FFT on regular optical networks[END_REF][START_REF] Chen | Wavelength assignment for directional hypercube communications on a class of WDM optical networks[END_REF][START_REF] Libeskind-Hadas | Multicast routing and wavelength assignment in multihop optical networks[END_REF][START_REF] Yu | Routing and wavelength assignment for 3-aryn-cube in array-based optical network[END_REF][START_REF] Yu | Routing and wavelength assignment for 3-ary n-cube communication patterns in linear array optical networks for n communication rounds[END_REF][START_REF] Yuan | Optimal routing and channel assignments for hypercube communication on optical mesh-like processor arrays[END_REF].

The exchanged hypercube is a link-diluted variation of the hypercube nework, proposed by Loh et al [START_REF] Loh | The exchanged hypercube[END_REF], with numerous desirable properties, such as lower diameter and better cost effectiveness. Some related works on exchange hypercubes, such as the domination number [START_REF] Klav_Zar | The domination number of exchanged hypercubes[END_REF], the connectivity [START_REF] Ma | The connectivity of exchanged hypercubes[END_REF], the super connectivity [START_REF] Ma | The super connectivity of exchanged hypercubes[END_REF], and faulttolerance measures [START_REF] Li | Generalized measures of fault tolerance in exchanged hypercubes[END_REF] have been investigated.

The rest of this paper is organized as follows. In Section 2, we introduce some preliminaries of exchanged hypercubes and the congestion of embedding schemes. In Section 3, a lower bound of the number of required wavelengths for realizing exchanged hypercubes communication patterns on linear arrays is obtained. In Section 4, we propose an embedding scheme and an optimal wavelength assignment algorithm. Finally, we conclude the paper in Section 5.

Preliminaries

In this section, we introduce some preliminaries of exchanged hypercubes and the congestion of embedding schemes.

The exchanged hypercube

Let n be a positive integer. The n-dimensional hypercube (or n-cube for short) Q n is the graph with vertex set {0, 1} Let k > 1 and u = u k-1 …u 0 ∈{0, 1} n be a binary string. We use u j:i to denote the substring u j u j-1 …u i of u for 0 ≤ i ≤ j < k.

Definition 2.1 ([9]

). The vertex set V of exchanged hypercube EH(s, t) is the set {u s+t u s+t-1 …u 0 u i ∈{0, 1} for 0 ≤ i ≤ s+t}. Let u = u s+t …u 0 and v = v s+t …v 0 be two vertices in EH(s, t). There is an edge (u, v) in EH(s, t) if and only if (u, v) is in one of the following sets:

E 1 = {(u, v)  u 0 ≠ v 0 , u i = v i for 0 ≤ i ≤ s+t }. E 2 = {(u, v)  u 0 = v 0 = 1, H(u, v) = 1 with u i ≠ v i for some 1 ≤ i ≤ t }, and E 3 = {(u, v)  u 0 = v 0 = 0, H(u, v) = 1 with u i ≠ v i for some t+1 ≤ i ≤ s+t }.
Let EH i (s, t) be the subgraph of EH(s, t) induced by the edges in E i for i ∈{1, 2, 3}. Clearly, EH(s, t) contains 2 s+t+1 nodes and is a spanning subgraph of hypercube Q s+t+1 . For u ∈ V(EH(s, t)), if u 0 = 0, then the degree of u is s+1; otherwise, the degree of u is t+1. Fig. 1 depicts EH(1, 2) which is a spanning subgraph of Q 4 . An edge with a label i for i ∈{1, 2, 3} is in edge set E i . We can see that each node u in EH(1, 2) with u 0 = 0 is of degree 2 and all the other nodes are of degree 3.

Lemma 2.2 ([9]). EH(s, t) is isomorphic to EH(t, s).

By Lemma 2.2, hereafter, we may assume without loss of generality that s ≤ t.

Proposition 2.3 ([6]

). EH 2 (s, t) (respectively, EH 3 (s, t)) contains 2 s (respectively, 2 t ) copies of Q t (respectively, Q s ) in which any two distinct copies of Q t (respectively, Q s ) are disjoint. Moreover, EH 1 (s, t) forms a perfect matching between nodes in EH 2 (s, t) and EH 3 (s, t).

Denote by

1 : + + t t s u t Q for the Q t in EH 2 (s, t) in which all vertices u ∈ Q t have the same bits in u s+t:t+1 . Similarly, 1 : t u s Q denotes those Q s in EH 3 (s, t) for all vertices u ∈ Q s hav-
ing the same bits in u t:1 . For brevity, Fig. 2. The two subgraphs of EH(1, 2).

The congestion

Let G = (V 1 , E 1 )
Φ = Φ ∈ The congestion of G in H under Φ, denoted by c g (G, H), is defined as: ) , , ( min ) , , ( Φ = Φ
Φ H G c H G c p g
Let λ(G, H) stand for the number of required wavelengths for realizing communication pattern G on WMN optical network H by embedding scheme Φ. Lemma 2.5 shows that c g (G, H) is a lower bound of λ Φ (G, H). Lemma 2.5 ([1, 3, 12]). λ Φ (G, H) ≥ c g (G, H).

In this paper, we consider that the guest graph is EH(s, t) and the host graph is a linear array L n , where n = s+t+1 and L n is a path of 2 n nodes. We label the nodes (respectively, the edges) in L n from 1 to 2 n (respectively, from e 1 to 1 2n e ) in consecutive order. For example, a linear array L 3 is shown in Fig. 3. Given an embedding scheme, each node u in V(EH(s, t)) will be assigned a distinct number in {1,…,2 n }. The node assigned number i is then embedded to the node i in L n . 3 A lower bound of λ Φ (EH(s,t), L n )

If u = u s+t …u t+1 u t …u 1 u 0 is a node in Q i t (0 ≤ i ≤ 2 s-1
) and the decimal value of u t:1 is j, then we also use q i,j t to denote node u. Let R i t stand for the subgraph of EH(s, t) induced by the nodes in Q i t and all nodes in EH 1 adjacent to some vertex in Q i t .

Lemma

+ ≥ Φ - + Φ ≤ ≤ + . Theorem 3.4 c g (EH(s, t), L n ) ≥ 2 s+t-1 + 2 t /3.

Lemma 3.5

The number of required wavelengths to realize EH(s, t) communication patterns on linear array L n is not less than 2 s+t-1 + 2 t /3. Proof. By Lemma 2.5 and Theorem 3.4, the lemma is thus proved.

Q.E.D.

Optimal wavelength assignment for realizing EH(s, t) on L n

In this section, we first derive an embedding scheme, and then describe a routing and wavelength assignment algorithm. Let u = u s+t …u t+1 u t …u 1 u 0 be a node in V(EH(s, t)). We partition V(EH(s, t)) into eight disjoint subsets as follows: Clearly, the subgraph induced by S i (1 ≤ i ≤ 4) comprises 2 s-1 disjoint (t-1)-cubes, and the subgraph induced by S i (5 ≤ i ≤ 8) comprises 2 t-1 disjoint (s-1)-cubes. If s > 2 for the subgraph induced by S m (1 ≤ m ≤ 4), we denote the (t-1)-cube by Q m,i t-1 where i (0 ≤ i ≤ 2 s-1 -1) is the decimal value of u s+t:t+2 , and the node u in Q m,i t-1 is represented by q m,i,j t-1 , where j (0 ≤ j ≤ 2 t-1 -1) is the decimal value of u t:2 . Otherwise, if s = 1, the (t-1)-cubeis denoted by Q m,0 t-1 , and the node u in Q m,0 t-1 is denoted by q m,0,j t-1 , where j (0

≤ i ≤ 2 t-1 -1) is the decimal value of u t:2 .
Similarly, if s > 2, for the subgraph induced by S m (5 ≤ m ≤ 8), we denote the (s-1)cube by Q m,i s-1 , where i (0 ≤i ≤2 t-1 -1) is the decimal value of u t:2 , and the node in

Q m,i s-1
is represented by q m,i,j s-1 , where j (0 ≤ j ≤ 2 s-1 -1) is the decimal value of u s+t:t+2 . Otherwise if s = 1, the 0-cube with decimal value i (0 ≤ i ≤ 2 t-1 -1) in substring u t:2 is denoted as Q m,i 0 , and the only node in Q m,i 0 is enoted as q m,i,j 0 . Let u i = u i s+t …u i t+1 u i t …u i 1 u i 0 be a node in S i (1 ≤ i ≤ 8), and let v = v s+t-3 …v t+1 v t … v 1 v 0 be a binary string of length s + t -2 with decimal value x (0 ≤ x ≤ 2 s+t-2 -1). If s > 2, then let u i s+t:t+2 = v s+t-3:t-1 and u i t:2 = v t-2:0 ; otherwise, for the case s = 1, let u i t:2 = v t-2:0 . We can find that the nodes u i (1 ≤ i ≤ 8) form two reversed direction cycles, denoted by cycle 1 (x) and cycle 2 (x), respectively. Fig. 4 shows the two reversed direction cycles.

Fig. 4. Two reversed direction cycles in EH(s, t).

An embedding scheme α, which assign numbers to the nodes in EH(s, t), is shown in Table 1.

Table 1. An embedding scheme α.

Embedding scheme α

Input: An exchange hypercube EH(s, t). Output: The assigned number NUM(u), u ∈ V(EH(s, t)). begin

Step 1. Set k = 1;

Step 2. For each node u ∈ EH(s, t), set NUM(u) = NULL;

Step 3. For m = 1 to 4 For i = 0 to 2 s-1 -1

For j = 0 to 2 t-1 -1 NUM(q m,i,j t-1 ) = k; k = k + 1;

Step 4. For m = 5 to 8 For i = 0 to 2 t-1 -1

For j = 0 to 2 s-1 -1 NUM(q m,i,j s-1 ) = k; k = k + 1; end Property 4.1 In the embedding scheme α, if x ≠ y or m 1 ≠ m 2 or i 1 ≠ i 2 , then the nodes in Q m1,i1

x-1 and the nodes in Q m2,i2 y-1 are embedded into two disjoint linear subarrays of L n .

Proof. This property is clear from the embedding scheme α.

Q.E.D.

Fig. 5 shows the numbers assigned to the nodes in EH(1, 2) by the embedding scheme α. A routing and wavelength assignment algorithm β for realizing EH(s, t) communication patterns on L n is shown Table 2. Table 2. A routing and wavelength assignment algorithm β.

Routing and wavelength assignment algorithm β

Input: An exchange hypercube EH(s, t), and the assigned number NUM(u), u ∈ V(EH(s, t)). Output: The assigned number NUM(u), u ∈ V(EH(s, t)). begin Step 1. For x = 0 to 2 s+t-2 -1 assign 1 unused wavelength to link e in cycle 1 (x); assign 1 unused wavelength to link e in cycle 2 (x); Step 2. For m = 1 to 4

For i = 0 to 2 s-1 -1 Call Algorithm 1 in [START_REF] Chen | Routing and wavelength assignment for hypercube in array-based WDM optical networks[END_REF] to assign wavelengths to links in Q m,i t-1 .

Step 2. For m = 5 to 8

For i = 0 to 2 t-1 -1 Call Algorithm 1 in [START_REF] Chen | Routing and wavelength assignment for hypercube in array-based WDM optical networks[END_REF] to assign wavelengths to links in Q m,i s-1 . end Theorem 4.2 The optimal number of required wavelengths to realize EH(s, t) communication patterns on L n is 2 s+t-1 + 2 t /3. Proof. It is clear that Algorithm β considers all links in EH(s, t). In Step 1, we have that links on 2 s+t-1 directed cycles are assigned wavelengths, and links on each directed cycle are assigned 1 unused wavelength. Hence, 2 s+t-1 wavelengths are assigned in this step. In Step 2 (respectively, Step 3), Algorithm 1 is invoked to assign wavelengths to links in Q m,i t-1 (respectively, Q m,i s-1 ). According to the results in [START_REF] Chen | Routing and wavelength assignment for hypercube in array-based WDM optical networks[END_REF], it follows that 2 t /3 (respectively 2 s /3) wavelengths are required for each Q m,i t-1 (respectively, Q m,i s-1 ) in Step 2 (respectively, Step 3). By Property 4.1, the wavelengths assigned to links in each Q m,i s-1 and Q m,i t-1 can be reused, and hence, Steps 2 and 3 require 2 t /3 wavelengths. It is obvious that Algorithm β requires 2 s+t-1 + 2 t /3 wavelengths. By Lemma 3.5, an optimal wavelength assignment is achieved. This completes the proof.

Q.E.D.

Fig. 6 shows the wavelengths assigned to the links in EH(1, 2) by the routing and wavelength assignment algorithm β.

Concluding remarks

In this paper, we study the optimal wavelength assignment for realizing the exchanged hypercube EH(s, t) communication patterns on linear array WDM optical network L n by proving that c g (EH(s, t), L n ) > 2 s+t-1 +2 t /3. We also design an embedding scheme and a routing and wavelength assignment algorithm which assigns the optimal number of wavelengths.

For the case when s = t, the exchanged hypercube is reduce to the dual-cubes [START_REF] Li | Efficient collective communications in dual-cube[END_REF], and 2 2s-1 + 2 s /3 wavelengths are required when the guest network is the dual-cube. For a future research, it is worthwhile to parallelize the wavelength assignment algorithm proposed in this paper, and consider the routing and wavelengths assignment issues for other types of communication patterns, such as, crossed cubes, twisted cubes, recursive circulants, etc.

Fig. 1 .

 1 Fig. 1. An exchange hypercube EH(1, 2).

Q

  are also denoted by Q x t and Q y s , respectively, where x and y are the decimal values of u s+t:t+1 and u t:1 , respectively. Fig.2(a) and (b) show the two subgraphs EH 2 (1,2) and EH 3 (1, 2), respectively. Note that EH 2 (1, 2) contains Q 0 2 and Q 1 2 while EH 3 (1, 2) contains Q 0 1 , Q 1 1 , Q2 1 , and Q 3 1 .

Fig. 3 .

 3 Fig. 3. A linear array L 3 .

S 1 =

 1 {u : u t+1 = 0; u 1 = 0 and u 0 = 1}, S 2 = {u : u t+1 = 0; u 1 = 1 and u 0 = 1}, S 3 = {u : u t+1 = 1; u 1 = 0 and u 0 = 1}, S 4 = {u : u t+1 = 1; u 1 = 1 and u 0 = 1}, S 5 = {u : u t+1 = 0; u 1 = 0 and u 0 = 0}, S 6 = {u : u t+1 = 1; u 1 = 0 and u 0 = 0}, S 7 = {u : u t+1 = 0; u 1 = 1 and u 0 = 0}, and S 8 = {u : u t+1 = 1; u 1 = 1 and u 0 = 0}.

Fig. 5 .

 5 Fig. 5. The numbers assigned to nodes in EH(1, 2).

Fig. 6 .

 6 Fig. 6. The wavelengths assigned to links in EH(1, 2).

  Two vertices (strings) u and v in Q n are adjacent if and only if they differ in exactly one coordinate. Let H(u, v) denote the Hamming distance between u and v, namely the number of coordinates in which u and v are different. Thus two vertices u and v in Q n are adjacent if and only if H(u, v) = 1.

n .

  be the guest graph and H = (V 2 , E 2 ) the host graph, where V 1  = V 2 . An embedding scheme of G in H is an ordered pair Φ = (Ψ, Ω), where Ψ is a bijection from V 1 to V2, Ω is a mapping from E 1 to a set of paths in H such that, for every edge e = (u, v) ∈ E 1 , there is a path Ω(e) from Ψ(u) to Ψ(v) in H. The congestion of an edge e ∈ E 2 under embedding scheme Φ of G in H, denoted by c e (G, H, Φ, e), is the number of paths Ω(e ' ) for all e' ∈ E 1 passing through e, namely,

	Definition 2.4. ) , , , ( max ) , , ( 2 e H G c H G c e E e p

c e (G, H, Φ, e) = {e ' : e ∈ E(Ω(e ' )), e' ∈ E 1 and e ∈ E 2 }

The congestion of G in H under Φ, denoted by c p (G, H, Φ), is defined as:

  3.1The congestion of embedding the nodes in Q i t to a linear subarray L t of L n is 2 t + 2 t /3. If Φ is an optimal embedding scheme, then Φ will embed nodes in a subgraph of EH(s, t), which is isomorphic to R i t , to nodes 1, 2,…2 t+1 of L n .

	Property 3.2 Lemma 3.3 e i	max 2 1 ; i t	1	{ min c e	(	EH	(	s	,	t	),	n L	,	,	i e	)}	2	s	t	1	  3 / 2 t