Sylvain Jay
email: sylvain.jay@fresnel.fr

Member, IEEE Mireille Olwen Guillaume
email: mireille.guillaume@fresnel.fr

Jacques Blanc Talon

Member, IEEE Jacques Blanc-Talon
email: mris.blanc-talon@sfr.fr

Underwater Target Detection With Hyperspectral Data: Solutions for Both Known and Unknown Water Quality

Keywords: Hyperspectral remote sensing, maximum likelihood estimation, underwater object detection

Underwater target detection with hyperspectral data : solutions for both known and unknown water quality",

I. INTRODUCTION

M APPING coastal water quality or benthic cover from remote-sensed data is a very challenging and widely studied problem: many authors have proved hyperspectral imagery could enable to retrieve efficiently bathymetry, bottom spectra or water constituents [START_REF] Lee | Hyperspectral remote sensing for shallow waters. II. deriving bottom depths and water properties by optimization[END_REF]- [START_REF] Vahtmäe | Feasibility of hyperspectral remote sensing for mapping benthic macroalgal cover in turbid coastal waters -a baltic sea case study[END_REF]. However, the detection problem in such a context remains open. For example, in [START_REF] Rodríguez-Díaz | Subsurface detection of coral reefs in shallow waters using hyperspectral data[END_REF], the purpose of the authors is more to classify pixels rather to detect target pixels among the background. In [START_REF] Masalmah | Clutter modeling for subsurface detection in hyperspectral imagery using markov random fields[END_REF], they consider only very shallow waters and do not take into account the attenuation by the water layer. However, underwater target detection from remote-sensing techniques is a very important problem; whether on a military or a civilian aspect, there are many challenges. In hyperspectral imagery, many detection algorithms have been developed during the last decades to meet many practical situations. Thereby, according to the starting hypotheses, each filter has several specificities: for example, when the mean spectral vector and covariance matrix of the background are available, the Matched Filter (MF) is developed, whereas the Adaptive Matched Filter (AMF) [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF] and Kelly Generalized Likelihood Ratio Test (GLRT) [START_REF] Kelly | An adaptive detection algorithm[END_REF] are convenient when these parameters have to be estimated. The Adaptive Cosine/Coherence Estimator (ACE) is very attractive in terms of false alarm stability because of its scaling invariance properties between training and test data [START_REF] Kraut | Adaptive subspace detectors[END_REF]. Unfortunately, the efficiency of these supervised detectors decreases in the case of underwater target. Actually, bottom and target spectra are considerably warped by the water column, particularly in the visible and the near-infrared range. This distortion depends not only on the depth but also on the water quality [START_REF] Maritorena | Diffuse reflectance of oceanic shallow waters -influence of water depth and bottom albedo[END_REF]- [START_REF] Brando | A physics based retrieval and quality assessment of bathymetry from suboptimal hyperspectral data[END_REF]. In order to keep good detection performances for such scenes, the idea is to insert a bathymetric model in the starting hypotheses before developing the filters, to correct these spectral distortions.

In section II, we present the chosen bathymetric model of reflectance, before using it in section III to develop new versions of three of the most popular filters, namely MF, AMF and ACE. In this section, we also propose a GLRT-based filter, for the case of unknown depth and concentrations of water constituents. This filter is designed to be self-sufficient for underwater target detection (e.g. it does not need any knowledge on water characteristics). We give the results of this estimation method in section IV, and we compare the new bathymetric detectors with usual filters applied after having performed an inversion process. Both simulated and real data are tested in this paper.

II. BATHYMETRIC MODEL OF REFLECTANCE A. General equation

As previously mentioned, many authors have focused on the ability of hyperspectral data to evaluate water characteristics. Therefore, many models have been proposed for subsurface reflectance, depending for example on the aim of the study or on the zone to consider: in order to estimate phytoplankton concentration in open sea or in such turbid waters that bottom has no influence on the subsurface reflectance, a first-order model is proposed in [START_REF] Kutser | A hyperspectral model for interpretation of passive optical remote sensing data from turbid lakes[END_REF]. If the aim is to retrieve bathymetry, the role of the bottom must be considered and more complex models are developed as in [START_REF] Maritorena | Diffuse reflectance of oceanic shallow waters -influence of water depth and bottom albedo[END_REF]. Firstly, since we study underwater target detection in optically shallow waters, we have chosen the well-known model of Maritorena et al. [START_REF] Maritorena | Diffuse reflectance of oceanic shallow waters -influence of water depth and bottom albedo[END_REF], transformed in [START_REF] Lee | Hyperspectral remote sensing for shallow waters. I. a semianalytical model[END_REF] to the subsurface remote-sensing reflectance as follows: where r ∞ (λ) represents the subsurface remote-sensing reflectance over an optically deep water column, r B (λ) the bottom albedo, k(λ) the attenuation coefficient and H the depth (each parameter and acronym is defined in Table I). The bottom is assumed to be a Lambertian reflector. The first term in (1) is the water column contribution (predominant in deep waters) whereas the second one is the bottom contribution.

r(λ) = r ∞ (λ) ( 1 -e -2k(λ)H ) + r B (λ) π e -2k(λ)H (1)
A more realistic model is presented in [START_REF] Lee | Hyperspectral remote sensing for shallow waters. II. deriving bottom depths and water properties by optimization[END_REF] and has been widely used during the last decade [START_REF] Brando | A physics based retrieval and quality assessment of bathymetry from suboptimal hyperspectral data[END_REF], [START_REF] Xiu | Deriving depths of deep chlorophyll maximum and water inherent optical properties: a regional model[END_REF]- [START_REF] Volpe | Remote sensing retrieval of suspended sediment concentration in shallow waters[END_REF]; in this expression, the attenuation coefficient k(λ) is different in the upwelling and downwelling directions (resp. k u (λ) and k d (λ)). One also separates the photons that never interact with the bottom from those that interact with the bottom:

r(λ) = r ∞ (λ) ( 1 -e -(k d (λ)+k c u (λ))H ) + r B (λ) π e -(k d (λ)+k b u (λ))H . ( 2
)
In this article, the simplified model (denoted as M1) is used for simulations; the second one (denoted as M2) is used with real data in order to assess the impact of a more realistic model on estimation and detection. The theoretical developments are made with M1, but can easily be written with M2.

For convenience, we ignore the effects of the water-air interface, which can be represented by the empirical relation of Lee et al. [START_REF] Lee | Hyperspectral remote sensing for shallow waters. I. a semianalytical model[END_REF] and connected after the detection process, as well as the atmosphere effects [START_REF] Kruse | Comparison of ATREM, ACORN, and FLAASH atmospheric corrections using low-altitude aviris data of boulder, co[END_REF].

B. Modeling of physical parameters

The water quality plays an important role in the light attenuation through two physical parameters, respectively the absorption and the backscattering coefficients a(λ) and b b (λ). They are mainly dependent on three optically active constituents, phytoplankton pigments, colored dissolved organic matter (or gelbstoff) and non-algal particles. Usually, the absorption and backscattering coefficients are described as the sum of the contributions of each component. In [START_REF] Brando | A physics based retrieval and quality assessment of bathymetry from suboptimal hyperspectral data[END_REF], Brando et al. express them as following (neglecting the CDOM influence on backscattering):

a(λ) = a w (λ) + C φ a * φ (λ) + C CDOM e -SCDOM (λ-λ0) + C N AP a * N AP (λ 0 )e -SNAP (λ-λ0) (3) b b (λ) = b b,w (λ) + C φ b * b,φ (λ 1 ) ( λ 1 λ ) Yφ + C N AP b * b,N AP (λ 1 ) ( λ 1 λ ) Y N AP (4)
where the index w is related to pure water, φ to phytoplankton, CDOM to gelbstoff and NAP to non-algal particles. a * j (λ) and b * b,j (λ) spectra are the specific absorption and backscattering coefficients. Reference wavelengths are fixed: λ 0 = 440 nm and λ 1 = 542 nm. We get the a w (λ), a * φ (λ) and b b,w (λ) spectra from [START_REF] Segelstein | The complex refractive index of pure water[END_REF], [START_REF] Bricaud | Variability in the chlorophyll-specific absorption coefficients of natural phytoplankton: analysis and parameterization[END_REF] and [START_REF] Richardson | Remote Sensing of Aquatic Coastal Ecosystem Processes : Science and Management Applications[END_REF] respectively. The remaining parameters can be found in [START_REF] Brando | A physics based retrieval and quality assessment of bathymetry from suboptimal hyperspectral data[END_REF]. The absorption and backscattering coefficients affect the subsurface reflectance in (1) through the reflectance of an infinitely deep water column r rs,∞ (λ) and, of course, the attenuation coefficients. There are several relations describing these links; for r rs,∞ (λ), they can be linear or quadratic with respect to [START_REF] Lee | Hyperspectral remote sensing for shallow waters. II. deriving bottom depths and water properties by optimization[END_REF], [START_REF] Kutser | A hyperspectral model for interpretation of passive optical remote sensing data from turbid lakes[END_REF], or even more complex as in [START_REF] Albert | An analytical model for subsurface irradiance and remote sensing reflectance in deep and shallow case-2 waters[END_REF], where Albert and Mobley take into account the surface wind, the viewing angle and the solar zenith angle. We choose expressions of intermediate complexity for both parameters: k(λ) is taken from [START_REF] Legleiter | A forward image model for passive optical remote sensing of river bathymetry[END_REF], whereas r rs,∞ (λ), k d (λ), k b u (λ) and k c u (λ) are taken from [START_REF] Lee | Hyperspectral remote sensing for shallow waters. II. deriving bottom depths and water properties by optimization[END_REF].

b b (λ) b b (λ)+a(λ)

C. Algebraic expression of reflectance vector

In order to use [START_REF] Lee | Hyperspectral remote sensing for shallow waters. II. deriving bottom depths and water properties by optimization[END_REF] in the context of hyperspectral imagery, we use a vector representation on the set {λ 1 , λ 2 , ..., λ L }, to lead to the following matrix form:

ρ = Kρ 0 (5) with ρ = r -r ∞ , ρ 0 = 1 π r B -r ∞ and K = diag(e -2k(λ1)H , ..., e -2k(λL)H ).
In [START_REF] Rodríguez-Díaz | Subsurface detection of coral reefs in shallow waters using hyperspectral data[END_REF], we can see clearly that the vector ρ just below the water surface is the vector ρ 0 calculated just above the bottom and attenuated by the matrix K. The latter describes the attenuation and contains information on depth H and inherent water properties a(λ) and b b (λ).

III. DEVELOPMENT OF BATHYMETRIC DETECTORS

In this article, we develop target detection algorithms using a statistical approach. Even though the assumption of normality is not always respected by hyperspectral data, many studies for multivariate Gaussian distributions have been led with success. So, we assume the spectral variability of ρ is Gaussian and models only the intrinsic variability of the class of the material. Actually, we do not consider mixed pixels but only pure ones. However, the bathymetric model ( 5) is linear and therefore, the linear mixing of bottom pixels is not affected by the water attenuation, and the abundances are preserved; the subpixel detection algorithms presented in [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF] could be then developed in this context. Two sets of data are used: training pixels are supposed "targetfree" whereas a test pixel can represent either background or target. To decide whether the test pixel ρ belongs to the target class (hypothesis H 1 ) or the background class (hypothesis H 0 ), we derive a contrast function D from the likelihood ratio given by the ratio of the two conditional probability density functions under each hypothesis:

D(ρ) = P (ρ|H 1 ) P (ρ|H 0 ) H1 ≷ H0 µ. ( 6
)
If D(ρ) is larger than µ, ρ is considered as a target pixel.

The likelihood ratio can lead to various functions according to the corresponding hypotheses and a priori knowledge. In the following sections, we transpose the well-known MF, AMF and ACE detectors in the context of bathymetry. We could mention the famous Kelly GLRT [START_REF] Kelly | An adaptive detection algorithm[END_REF], but it turns out that, if there are enough training pixels, the AMF and the Kelly GLRT give some very similar results [START_REF] Kraut | Adaptive subspace detectors[END_REF], [START_REF] Robey | A CFAR adaptive matched filter detector[END_REF]. Thus, we choose to work with a number of training pixels at least five times greater than the number of bands, in order to be sure the performances of both filters are equivalent. However, it's helpful to keep in mind that the Kelly GLRT is more efficient than the AMF for low training data support. Furthermore, we propose a new GLRT-based detector if some water characteristics are unknown, in order to develop underwater target detection without any a priori knowledge on the water column.

A. Bathymetric matched filter (BMF)

For each test pixel, we consider the two following hypotheses:

• H 0 : background pixel ⇒ ρ ∼ N (µ b , Γ b ) • H 1 : target pixel ⇒ ρ ∼ N (µ t , Γ t ) Note that µ b = K (r b -r ∞ ) and µ t = K(r t -r ∞ )
where r b and r t are the bottom and target spectra respectively. Each parameter is assumed to be known. The logarithm of the likelihood ratio (6) written with Gaussian conditional probability density functions leads to the following linear bathymetric matched filter (assuming that the background and target classes share the same covariance matrix Γ [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF]):

d BM F (ρ) = (µ t -µ b ) t Γ -1 (ρ -µ b ). (7)

B. Bathymetric adaptive matched filter (BAMF)

The previous BMF assumes all parameters are known, including covariance matrix. In practical situations, the latter must be estimated from data. An adaptive algorithm has been developed in the original paper of Robey et al. [START_REF] Robey | A CFAR adaptive matched filter detector[END_REF], and transposed in the hyperspectral context in [START_REF] Farrell | Analysis of modeling, training, and dimension reduction -approaches for target detection in hyperspectral imagery[END_REF]. The hypotheses differ from the BMF in the sense that, under hypothesis H 0 , the centered spectral vectors x = ρµ b are distributed as a Gaussian random vector with mean 0 and covariance matrix Γ, and under hypothesis H 1 , it is assumed to have mean b(µ t -µ b ) and covariance matrix Γ. It means that the difference µ t -µ b can be known up to a multiplicative factor b. Therefore, we have:

• H 0 : background pixel ⇒ x ∼ N (0, Γ) • H 1 : target pixel ⇒ x ∼ N (b(µ t -µ b ), Γ)
The procedure, as described in [START_REF] Robey | A CFAR adaptive matched filter detector[END_REF], leads to the following expression for the bathymetric adaptive matched filter):

d BAM F (ρ) = [ (µ t -µ b ) t Γ -1 (ρ -µ b ) ] 2 (µ t -µ b ) t Γ -1 (µ t -µ b ) . ( 8
)

C. Bathymetric Adaptive Coherence/Cosine Estimator (BACE)

Another filter currently used is the Adaptive Coherence/Cosine Estimator. Although one of its qualities is to deal better than the above detectors with subpixel situations, it also has an excellent power of discrimination between target and other spectra present in the scene. In [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF], the authors express the ACE filter as the square cosine of the angle between the test pixel and the target pixel in the whitened coordinate space. The corresponding bathymetric filter is:

d BACE (ρ) = cos 2 ( Γ -1 2 (ρ -µ b ), Γ -1 2 (µ t -µ b ) ) (9)
d BACE (ρ) = [ (µ t -µ b ) t Γ -1 (ρ -µ b ) ] 2 * 1 [ (µ t -µ b ) t Γ -1 (µ t -µ b ) ] [ (ρ -µ b ) t Γ -1 (ρ -µ b )
] .

(

) 10
This detector has been also developed with a GLRT approach in [START_REF] Kraut | The cfar adaptive subspace detector is a scale-invariant glrt[END_REF] by Kraut and Scharf.

D. A solution in case of unknown water parameters

Usually, the physical parameters such as depth or phytoplankton concentration are unknown a priori. But they are needed in order to correct the target spectrum before detection (from r t to µ t ). Therefore, the objective is to develop a selfsufficient filter. We estimate the unknown set of parameters θ using a maximum likelihood approach and developing the GLRT with all the data (test pixel ρ and training pixels χ = {ρ 1 , ..., ρ N }).

1) A GLRT-based bathymetric filter (GBF): For example, if both covariance matrix and depth are unknown, we consider the following GLRT:

D GBF (ρ) = max θ P (ρ, χ; θ|H 1 ) max θ P (ρ, χ; θ|H 0 ) H1 ≷ H0 µ,
with θ = {Γ, H}. [START_REF] Lee | Hyperspectral remote sensing for shallow waters. I. a semianalytical model[END_REF] The development of this detector is showed in the appendix, and leads to:

D GBF (ρ) = 1 + (ρ -µ b ( H)) t S( H) -1 (ρ -µ b ( H)) 1 + (ρ -µ t ( H)) t S( H) -1 (ρ -µ t ( H)) . ( 12
)
In the same way, we can also develop this GLRT-based bathymetric filter for depth and concentrations:

D GBF (ρ) = max θ P (ρ, χ; θ|H 1 ) max θ P (ρ, χ; θ|H 0 ) H1 ≷ H0 µ, with θ = {Γ, H, C φ , C CDOM , C N AP }. ( 13
)
2) Parameter estimation: In the appendix, we see that, for a large number of training pixels, estimating H under the hypotheses H 0 or H 1 is equivalent to estimate H only on training data. We can generalize for the estimation of θ ′ = {H, C φ , C CDOM , C N AP }. Therefore, in the case of unknown covariance, depth and concentrations, we obtain:

θ ′ = argmax θ ′ |S(θ ′ )| -1 ( 14
)
where

S(θ ′ ) = N ∑ i=1 (ρ i -µ b (θ ′ ))(ρ i -µ b (θ ′ )) t .
In this article, we optimize the likelihood function with a simple optimization method, the relaxation method that alternatively maximizes each parameter in order to find a maximum (we suppose that this function is at least locally convex around the maximum).

IV. EXPERIMENTS AND RESULTS

A. Data 1) Simulated data: In simulated images, the bottom is built from a linear mixture of three of the main components of sand, namely quartz, feldspar and mica. The spectra are taken from the USGS spectral library [START_REF] Clark | Usgs digital spectral library splib06a: U.s. geological survey, data series 231[END_REF]. The intra-class variability is modeled with a white Gaussian noise (σ = 0.02). Once attenuated by the water column, the same white Gaussian noise is added to model other sources of variability, including sensor noise. It may be interesting to introduce a Poisson noise instead in order to cope better with recent sensors, which attain the Poisson noise limit; however, we also want to describe the natural spectral variability that could be due to inhomogeneities of the water quality and that can be well modeled by a Gaussian distribution. Therefore, this modification may not improve the results. In Fig. 1, we show the simulated images and the target and bottom spectra used in this article over and under a 3m-high turbid water column. In the entire article, the water constituent concentrations of turbid water are taken from [START_REF] Brando | A physics based retrieval and quality assessment of bathymetry from suboptimal hyperspectral data[END_REF]. The wavelength range is fixed at 400-700 nm, with a spectral resolution of 5 nm (61 spectral bands). Herein, the training set contains 441 pixels in order to be sure that the BAMF and Kelly GLRT are equivalent, and that the ML estimates for physical parameters under hypotheses H 0 and H 1 are almost equal. The signal-to-noise ratio (SNR) in decibels is calculated as: SNR = 10 log

( ∑ i ρ 2 i / ∑ i n 2 i ) in which n i
is the simulated noise in pixel ρ i . The SNR varies with the additive sensor noise. 2) Real data: Real data have been collected by ACTIMAR 1as part of the HYPLITT project, supported by the Office for Advanced Research and Innovation (DGA/MRIS), France. The study site is located in Quiberon Peninsula, on the West coast of France (see Fig. 2). The spatial resolution is 0.4 m whereas the spectral resolution is about 4.5 nm. The ATCOR model [START_REF] Richter | Geo-atmospheric processing of airborne imaging spectrometry data, part.2: atmospheric/topographic correction[END_REF] is used to perform hyperspectral atmospheric correction. Finally, we use the model of Lee et al. [START_REF] Lee | Hyperspectral remote sensing for shallow waters. II. deriving bottom depths and water properties by optimization[END_REF] to describe the water-air interface crossing. In order to assess the detection performances, black and white tarpaulins have been placed on a sandy bottom at different depths: the depths in I1 and I2 are respectively 4.70 m and 6.70 m. The target size is 1 m 2 . Moreover, a greater range of depths (between 0.20 m and 16.80 m) is available in other images but unfortunately, the water quality has not been measured.

B. Comparison with classical methods

We compare our estimation method against the least square method, in which the error function is optimized by the Levenberg-Marquardt (LEVMAR) algorithm, as proposed in [START_REF] Albert | An analytical model for subsurface irradiance and remote sensing reflectance in deep and shallow case-2 waters[END_REF] and [START_REF] Xiu | Deriving depths of deep chlorophyll maximum and water inherent optical properties: a regional model[END_REF]. In order to assess the bathymetric filters, we firstly perform an inversion with the parameters estimated with the LEVMAR method, before applying the usual detectors (AMF, ACE) to the transformed images. The Matlab implementation of the LEVMAR algorithm is used.

C. Estimation results

The above-mentioned GLRT-based filter estimates the missing physical parameters before detection. Consequently, the The estimation results depend on the a priori knowledge on model parameters. So the precision is different if we estimate the covariance matrix and depth only, or covariance matrix and all the other physical parameters. In simulations, we estimate the covariance matrix, depth and concentrations at the same time, and therefore, we maximize the function ( 14) with the optimization scheme described in section III-D2. As a consequence, the errors would be lower if we had more a priori knowledge on model parameters. The results are obtained with training pixels χ = {ρ 1 , ..., ρ N } following the hypothesis H 0 of section III. In real scenarios, it seems feasible to obtain this training set since it comes to consider an area of 64 m 2 around the central pixel (with the spatial resolution of HYPLITT images). Near the sea shores, the depth do not vary very much within such small zones. Otherwise, it is still possible to reduce the size of the training set without lowering the results dramatically. In our future

Y = E χ [ ( Y -Y ) 2
] . Excepted for very shallow water (0.1 m), in which the water constituents do not have much influence on subsurface reflectance, we note that, on simulated data, the deeper the For optically deep water (H = 30 m in these conditions of experiments), the bottom has almost no more influence on subsurface reflectance and we have: ρ 0 ⇐⇒ r r ∞ . Therefore, r hardly does not depend on H any more and the estimation of the depth becomes tricky. On the contrary, the precision of the estimation of C φ , C CDOM and C N AP does not decrease beyond 20-30 m, since the reflectance remains equal to r ∞ . We also show the estimation robustness for high noise levels in Table III for H = 14 m. From an overall point of view, we see that the estimation errors remain moderate. As in previous tables, we note that the most sensitive parameter is C φ whereas H remains well estimated, even for SNR = 1 dB. Furthermore, we test this method on several real images presented in the previous section. Depths vary from 0.20 m to 16.80 m. It is good to notice that, contrary to simulated data, sandy bottoms are not perfectly plane, and these small variations of depth in each image definitely reduce the precision. Unfortunately, water characteristic measurements could not be realized for most of the data. The dry sand, wet sand and target spectra have been acquired on the beach before the flight; coupled with the depth H, they constitute the ground truth. Therefore, in the following, we focus on the results of depth estimation. We consider two images for each depth: their size vary from 200 to 1000 pixels. In Table IV, we show the relative RMSE of the estimation process for both maximum likelihood estimation (MLE) and LEVMAR methods and both models (1) and ( 2). Firstly, we can note that MLE performs better than LEVMAR in most cases: for example, at a depth of 2.83 m and using the M1 model, the RMSE for MLE is 22.4% whereas the RMSE for LEVMAR algorithm is 40.5%. For this water quality and the M1 model, depths lower than 12 m can be estimated with MLE with a RMSE lower than about 30%. In addition, using a more realistic model seems to be quite interesting for improving easily the estimation results, especially for MLE. The estimation errors are greater than in Table II, and depend on the locations. Several reasons can explain that: the bottom, depth or concentrations may change within the image. The sunglint and some modeling errors can also affect the precision. By way of example, we compare in Fig. 3 the simulated subsurface sand spectrum (obtained with the real sand spectrum corrected with the estimated parameters) with several observed spectra (H = 1.12 m) and we see that the simulated spectrum fits quite well the observed spectra. Therefore, the general trend emerging from these results indicates that this estimation process performs well: however, for very shallow waters, the effect of the water column is low and smothered by noise, and then it is difficult to estimate the water quality. For deep waters, the bottom does not have much influence on reflectance and the depth estimation is less precise. The resulting optimal range of depths is not fixed, and especially depends on the noise level, water quality and number of training pixels.

D. About the detection process

In real scenarios, the detection process could be as following: assuming that the target pixels are sparse enough and then, do not affect the estimation process, we estimate the parameters in the test pixel neighborhood, and calculate the output of the filter for each pixel of this neighborhood. In this section, we prove the interest of this underwater target detection method. In Fig. 4.a, we compare the performances of each filter for detecting galvanized metal (first-line targets in Fig. 1) under a 55m-high pure water column and a SNR of 5.6 dB. From an overall point of view, we can see clearly that the bathymetric filters outperform the LEVMAR-based methods. For a probability of false alarm (PFA) of 10 -3 , bathymetric filters detect around 80% of target pixels whereas LEVAMF or LEVACE do not detect any target pixels. These performances depend obviously on the differences between target and bottom spectra, and on the water quality. For more turbid water (see Fig. 4.b), the probability of detection is around 0.7 for a PFA of 10 -3 at 14 m and for SNR = 9.9 dB.

These previous curves show the performances for a fixed SNR. But it would be interesting to study the noise robustness of these detectors; indeed, we can fear that the GLRT-based filter is sensitive for low SNRs because the precision of the estimation decreases when the noise level increases. In Fig. 5, we plot the probability of detection of galvanized metal as a function of SNR for a given depth. We can see that BMF and GBF detect better than BAMF and BACE, even for noisy images. GBF is more robust than BAMF, as both filters are rather equivalent for high SNRs, whereas for low SNRs, its performances decrease slower than those of BAMF. Then, as previously showed in [START_REF] Jay | Underwater target detection with hyperspectral remote-sensing imagery[END_REF] with BMF, the parameter estimation does not have too much influence on detection performances. Besides, these results can be combined with those of Table III: we note that even if the RMSEs on C φ and C CDOM is respectively about 17% and 11% for SNR = 5 dB, GBF remains rather equivalent to the other filters. In such conditions, this GLRT-based filter is robust to estimation errors on physical parameters and can be used without any a priori knowledge on water column. We also implement our algorithms (built with the M2 model) on the real data presented in Fig. 2. We present the associated detection results for both I1 and I2 sites in Fig. 6. The global tendency confirms that the bathymetric filters outperform the LEVMAR-based detectors. For such depths, the white tarpaulins are well detected by each filter, as all target pixels are detected for a PFA of 10 -3 . Regarding the black tarpaulins that are more difficult to detect, BAMF and BACE perform slightly better than LEVAMF and LEVACE respectively in I1 and I2. In this case, we can also notice that GBF is the most efficient detector: for example, in I1 and I2 and for a PFA of 10 -3 , the PD is equal to 0.8 whereas it is lower than 0.2 for the other detectors. As mentioned previously, this confirms the fact that GBF overcomes the other filters in tricky conditions, despite the estimation errors on the depth and concentrations. Moreover, let note that in these images, the existing sunglint can affect the performance and therefore, a prior correction would certainly increase it. As far as the computation time is concerned, we compare the per-pixel processing time for both estimation and detection.

For LEVMAR-based methods and with a 3 Ghz processor, it is equal to 12 ms and 13 ms for I1 and I2 respectively, while our algorithm needs 4 ms and 2.4 ms. Of course, it is possible to improve easily our algorithm using a more efficient optimization scheme, such as gradient-based techniques, since the relaxation method is not very fast. Therefore, these last results confirm the interest of such bathymetric filters for the underwater target detection.

V. CONCLUSIONS

In this article, we develop a new approach for detecting underwater targets. We use a bathymetric model of subsurface reflectance, which is described as a weighted sum of the bottom and water column influences. This well-known physical model takes into account the water quality through the concentrations of three optically active constituents, e.g. the phytoplankton pigments, CDOM and non-algal particles. Then, using a statistical approach and considering a pixel as a random vector parameterized by the above-mentioned depth and concentrations, we develop new bathymetric detectors, based on usual ones. In real scenarios, most of the time, some of the water characteristics are unknown and therefore, we also propose the GBF, a GLRT-based filter that estimates these parameters on training pixels for a self-sufficient detection. This filter shows good robustness properties, and copes well with noisy scenarios and low contrasted targets. It is particularly adapted to difficult detection situations, such as underwater detection with sunglint effects. This estimation method is tested on both simulated and real data and assessed as a function of the depth and SNR. The global results are good, especially for depths between a few meters and about 20 m, and a SNR greater than 5 dB. The estimation of C φ is the most sensitive for usual values of depth and SNR, whereas the one of C N AP is not very sensitive. The estimation of H is quite good for such depths and SNRs, but is much more difficult for turbid waters deeper than 20 m, as the bottom has no more influence on subsurface reflectance. Besides, we observe the same phenomenon for the estimation of H on real data. On the other hand, the precision of concentration estimation does not decrease any more beyond 20 m since the reflectance remains equal the reflectance of an infinitely deep water column. However, we have to keep in mind that these depth values depends on the noise level, water quality and number of training pixels. Moreover, we see that the insertion of a more realistic bathy-metric model enables to increase easily the estimation results. The underwater detection algorithms are tested on simulated and real data. We see that the proposed detectors overcome the classical ones for deep areas, even when the parameter estimation is less accurate. On real hyperspectral images, good detection rates are obtained for depths between 4 m and 7 m. The misestimations of covariance matrix and physical parameters H, C φ , C CDOM and C N AP can impact the detection performances but still improve them with respect to usual filters, and therefore, allow the detection without any a priori knowledge on water column. In our future works, we will study the Cramer-Rao bounds for the estimation of model parameters. Under some hypotheses on the covariance structure, some works have already been done on simulated data [START_REF]Estimation of water column parameters with a maximum likelihood approach[END_REF]; we will then extend our methods to a more general case, and test them on real data. We will also introduce a parametric model for the depth in order to cope with varying depths in the studied area.

As Kelly does in [START_REF] Kelly | An adaptive detection algorithm[END_REF], we first estimate the covariance matrix under each hypothesis. For H 0 and H 1 , we obtain respectively:

Γ 0 = 1 N + 1 [ (ρ -µ b (H))(ρ -µ b (H)) t + S(H) ] (16
)

Γ 1 = 1 N + 1 [ (ρ -µ t (H))(ρ -µ y (H)) t + S(H) ] (17)
where S(H)

= N ∑ i=1 (ρ i -µ b (H))(ρ i -µ b (H)) t is obtained
from training data. The GLR becomes:

D GBF (ρ) = max H | Γ 1 | -1 max H | Γ 0 | -1 . ( 18
)
We can show that, if the number of training pixels N is large enough, optimizing the likelihood function under hypothesis H 0 is equivalent to optimize it under hypothesis H 1 or only on training data. This can be verified in Fig. 7: we plot the MSE of the difference between the depth estimate under hypothesis H 0 or H 1 , and the depth estimate H calculated only on training data, as a function of N. We can see that if N is about seven times greater than the number of bands, the MSE is lower than 10 -4 m 2 , e.g. the error is lower than 1 cm. Thus, we can consider the ML estimates H H1 and H H0 are almost equal to H, and with the lemma for the determinant calculation used by Kelly in [START_REF] Kelly | An adaptive detection algorithm[END_REF], we finally obtain the following GLR (for unknown covariance matrix and depth):

D GBF (ρ) = 1 + (ρ -µ b ( H)) t S( H) -1 (ρ -µ b ( H)) 1 + (ρ -µ t ( H)) t S( H) -1 (ρ -µ t ( H)) . (19)

Fig. 1 .

 1 Fig. 1. Composite images from hyperspectral simulated data and the corresponding observed reflectance spectra. Left column: bottom without water. Right column: underwater bottom (H = 3 m, Cφ = 0.7 µg.L -1 , C CDOM = 0.08 m -1 , C N AP = 2.8 mg.L -1 ).

Fig. 2 .

 2 Fig. 2. HYPLITT true colour composite image. Two black tarpaulins and two white tarpaulins were placed in I1 and I2. The depth is 4.70 m for I1 and 6.70 m for I2.

Fig. 3 .Fig. 4 .

 34 Fig.3. Comparison between the simulated and real subsurface sand spectra (H = 1.12 m, RMSE = 16.9%)

Fig. 5 .

 5 Fig. 5. Probability of detecting galvanized metal as a function of SNR (H = 14 m, Cφ = 0.7 µg.L -1 , C CDOM = 0.08 m -1 , C N AP = 2.8 mg.L -1 ). The PFA is fixed at 10 -4 .

Fig. 6 .

 6 Fig. 6. ROC curves describing the detection of: a) white tarpaulin in I1: H = 4.70 m, Cφ = 1.5 µg.L -1 , C CDOM = 0.09 m -1 , C N AP = 3 mg.L -1 ; b) white tarpaulin in I2: H = 6.70 m, Cφ = 2.4 µg.L -1 , C CDOM = 0.09 m -1 , C N AP < 2 mg.L -1 ; c) black tarpaulin in I1: H = 4.70 m, Cφ = 1.5 µg.L -1 , C CDOM = 0.09 m -1 , C N AP = 3 mg.L -1 ; d) black tarpaulin in I2: H = 6.70 m, Cφ = 2.4 µg.L -1 , C CDOM = 0.09 m -1 , C N AP < 2 mg.L -1 .

  APPENDIXIf some parameters (H, C φ , C CDOM or/and C N AP ) are unknown, we cannot use BMF, BAMF or BACE, because we need to know them for correcting the bottom and target spectra. For example, if both covariance matrix and depth are unknown, we have the following GLRT:D GBF (ρ) =max θ P (ρ, χ; θ|H 1 ) max θ P (ρ, χ; θ|H 0 ) H1 ≷ H0 µ, θ = {Γ, H}. (

Fig. 7 .

 7 Fig. 7. Maximum MSE of the difference between the depth estimate under hypothesis H 0 or H 1 , and the depth estimate H calculated only on training data, as a function of N (H = 10 m, turbid water).

TABLE II RESULTS

 II

		H = 0.1m H = 5m H = 10m H = 20m H = 30m
	H	1.18	0.33	0.50	1.62	26.81
	Cφ	32.53	2.95	6.66	18.76	19.02
	C CDOM	9.93	1.24	3.87	10.69	10.66
	C NAP	34.94	2.97	3.56	3.77	3.35
	works, we will try to model the bottom slope and insert it in
	the likelihood expression in order to get a better estimation of
	depth.

OF ML ESTIMATION IN THE CASE OF UNKNOWN COVARIANCE MATRIX, DEPTH AND CONCENTRATIONS (SIMULATED DATA): RELATIVE RMSE (%) AS A FUNCTION OF DEPTH (TURBID WATER) AND CALCULATED WITH 500 RUNS. The results are presented in Table

II

: we show the evolution of the relative Root Mean Square Error (RMSE) for the estimation of physical parameters, as a function of depth. The MSE for the parameter Y is given by: MSE

TABLE III RESULTS

 III OF ML ESTIMATION IN THE CASE OF UNKNOWN COVARIANCE MATRIX, DEPTH AND CONCENTRATIONS (SIMULATED DATA): RELATIVE RMSE (%) AS A FUNCTION OF SNR (H = 14 M) AND CALCULATED WITH 100 RUNS. the less precise the estimation of H, C φ and C CDOM is. Only the precision for C N AP seems to remain constant in this range of depths. The precision of the estimation of C φ and C CDOM goes down at H = 20 m (18% and 10% RMSE respectively) whereas H remains rather well estimated up to this depth (2% RMSE); besides, this is still true for very shallow water because, even if each other parameter are badly estimated, the RMSE for H is equal to 1.18%.

		SNR = 1dB SNR = 5dB SNR = 10dB SNR = 20dB
	H	3.00	1.18	0.76	0.35
	Cφ	51.65	17.44	10.77	5.81
	C CDOM	27.26	10.84	6.93	3.63
	C NAP	13.46	5.54	3.55	1.77
	water is,

TABLE IV RESULTS

 IV OF DEPTH ESTIMATION IN THE CASE OF UNKNOWN COVARIANCE MATRIX, DEPTH AND CONCENTRATIONS (REAL DATA): RELATIVE RMSE (%) FOR BOTH METHODS AND BOTH M1 AND M2 MODELS.

		H = 0.2m	H = 1.12m	H = 2.83m	H = 4.70m
		M1	M2	M1	M2	M1	M2	M1	M2
	MLE	22.3	7.4	25.6 16.9	22.4	9.0	30.8 15.6
	LEV	27.9	24.9 37.5 27.2	40.5 29.0 50.2 51.0
		H = 6.70m	H = 11.80m	H = 16.80m
		M1	M2	M1	M2	M1		M2
	MLE	31.3	26.7	10.3	17.8	528.1 122.7
	LEV	16.3	23.4 368.1	276.7	314.5 161.0

ACTIMAR is a firm specialised in a broad range of expertise areas within the marine environment sector, especially high resolution remote sensing applied to coastal zone management and deep ocean mapping.

This work was supported by the Office for Advanced Research and Innovation (DGA/MRIS).

S.


Jean-Max Redonnet

Walter Rubio

Gilles Dessein

Side milling of ruled surfaces: Optimum positioning of the milling cutter and calculation of interference

Keywords: five-axis machines, interferences, NURBS, ruled surfaces, side milling

Side milling is a process that enables machining time and thus costs, to be reduced. This type of machining is particularly well suited to ruled surfaces and all surfaces where one of the principal curvatures is very small compared to the tool radius and changes little oiler the entire surface. These surfaces must be treated with great care, as they are often on parts with high added value such as turbine blades, aircraft wings or helicoidal parts as used in fluid dynamics.

We then need to calculate and minimise interference that may arise if the I ruled surface cannot be developed. Whereas machining is usually carried out by setting the tool according to a rule, we suggest a new setting of the tool allowing interference to be reduced considerably. The computation algorithms for this setting were developed so as to be used in real-time by CAD/CAM software.

A comparative study of errors made with each type of setting is also presented. This shows a considerable reduction in errors when the setting we suggest is implemented.

Introduction

Side milling has been less frequently addressed by researchers than end milling (see [START_REF] Marciniak | Geometric Modelling for Numerically Controlled Machinig[END_REF][START_REF] Choi | Cutter location data optimization in 5-axis surface machining[END_REF][START_REF] Lee | Machined surface error analysis for 5-axis machining[END_REF][START_REF] Li | 5-axis machining of sculptured surfaces with a flat-end cutter[END_REF]), but it may be very useful, especially for ruled surfaces. We shall first briefly recapitulate the principle for ruled surfaces and settings that are currently adopted in existing CAD/CAM software before going on to introduce a new setting for the tool, allowing interference to be reduced considerably. The second part of the article is devoted to determining the parameters for this setting. We shall then make a comparative analysis based on the calculation of interference for each type of setting.

Ruled Surfaces

A ruled surface is generated by a set of straight lines based on two directrices.

The equation corresponding to the ruled surface in Fig. 1 is:

S(u, v) = (1 -v)C 0 (u) + vC 1 (u) (u, v) ∈ [0, 1] 2 C 1 (u) C 0 (u)
S(u, v)

P 0 P 1 (v = 0) (v = 1)
Figure 1: Ruled surface

"Standard" setting

With "standard" setting the tool axis is collinear with the rule under consideration. The tool can then 1 (a)

C 0 (u) ρ 0 C 1 (u)

ρ 1 P 0 N 1 N 0 α ε (b) C 0 (u) ρ 0 C 1 (u) ρ 1 N 1 P 0 N 0 N 2 α 0 α 1 ε ε Figure 2: "Standard" tool setting.
be placed at a tangent to one of the two directrices (Fig. 2a) or so as to distribute the error between its two extremities (solution proposed by Rubio [START_REF] Rubio | Calculation of Tool Paths for a Torus Mill on Free-Form Surfaces Five-Axis Machines with Detection and Elimination of Interference[END_REF][START_REF] Rubio | Génération de trajectoires du centre de l'outil pour l'usinage de surfaces complexes sur machines à trois et cinq axes[END_REF]; Fig. 2b).

According to [START_REF] Rubio | Calculation of Tool Paths for a Torus Mill on Free-Form Surfaces Five-Axis Machines with Detection and Elimination of Interference[END_REF], we then obtain in case (a):

ε = ρ 1 + R -R 2 + ρ 2 1 + 2Rρ 1 cos(α)
in case (b):

ε = ρ 0 + R -R 2 + ρ 2 0 + 2Rρ 0 cos(α 0 ) = ρ 1 + R -R 2 + ρ 2 1 + 2Rρ 1 cos(α 1 )
These equations can be used to determine the tool radius R needed to respect a given tolerance. In the case of infinite curvature radii (straight directrices), error calculation for the standard setting will give (see [START_REF] Rubio | Génération de trajectoires du centre de l'outil pour l'usinage de surfaces complexes sur machines à trois et cinq axes[END_REF]): in case (a):

ε = R(1 -cos(α)) in case (b): ε = R(1 -cos(α 0 ))

Optimum Setting

The setting we suggest is the fruit of a trade-off between two contradictory objectives:

1. To reduce interference by a more complex tool setting that therefore takes longer to calculate than the standard setting.

2. To be able to calculate this setting rapidly enough to enable real-time utilisation in CAD/CAM type software.

This optimum setting as described in Fig. 3 is based on 3 tangential points:

-Tangency of the tool lower generating line to the rule P 0 P 1 : point M 2 . -Tangency to the two directrices C 0 (u) and C 1 (u) on both sides of the extreme points of the rule: points M 0 and M 1 .

S(u, v) C 1 (u)
C 0 (u)

P 1 P 0 γ M 2 M 1 M 0 Figure 3: Optimum setting
Marciniak first proposed that the tool may be shifted along the rule (see [START_REF] Marciniak | Geometric Modelling for Numerically Controlled Machinig[END_REF] pp. 233,234), but the extra degree of freedom introduced in this way needs an additional drive surface to be eliminated. In [START_REF] Liu | Five-Axis NC cylindrical milling of sculptured surfaces[END_REF], Liu presents a method based on two offset points (the double point offset method) introducing an angle between the rule and the tool axis. In Section 4.2 we shall compare the results obtained using each method for the test surface presented in this paper.

With our setting, the tool position can thus be completely defined by two parameters: the angle γ made between the tool axis and the rule, and the position of the point M 2 on the rule P 0 P 1 (once the direction of the rotation axis is established).

Determining Setting Parameters 2.1 Preliminary Calculations

We first carry out the following preliminary calculations [START_REF] Carmo | Differential Geometry of Curves and Surfaces[END_REF][START_REF] Qiulin | Surface Engineering Geometry for Computer Aided Design and Manufacture[END_REF][START_REF] Faux | Computational Geometry for Design and Manufacture[END_REF]: The position of the rule is defined by u = u p . We then consider the vectors T 0 and T 1 tangent to the directrices C 0 (u) and C 1 (u) at the extremities of the rule P 0 and P 1 , and the vectors N 0 and N 1 perpendicular to the directrices C 0 (u) and C 1 (u) at P 0 and P 1 and also perpendicular to the rule P 0 P 1 , such that

N 0 = T 0 × P 0 P 1 and N 1 = T 1 × P 0 P 1
We establish R = tool radius h p = length of the rule at position

u p h 0 = distance P 0 , M 2 h 1 = distance P 1 , M 2 and obtain h 0 + h 1 = h p
We then define the following references (see Fig. 4):

1. Reference R 1 related to the rule:

-origin P 2 : at the centre of the rule -vector z 1 : along the direction of the rule -vector y 1 : along the direction of the bisector of the angle formed by the projections of N 0 and N 1 in the plane perpendicular to z 1 going through P 2 (vectorN 2 in Fig. 2b with α 0 = α 1 ) -vector x 1 : result of the vector product between y 1 and z 1

We can then calculate the matrices for passage between the overall reference in which the geometrical entities are defined and R 1 .

2. Reference R 2 related to the tool:

-origin P 2 : tangential point of the tool on the rule -vector y 2 : equal to y 1 -vector z 2 : projection of the tool axis on the plane P 2 , z 1 , x 1 -vector x 2 : result of the vector product between y 2 and z 2 In this reference, the equation for the tool will be:

a 0 L 0 h 0 θ 0 γ 0 M 2 P 0 x 2 z 2 x 1 z 1 x 2a z 2a plane of projection of C 0 (u) = plane of ellipse Ell 0 (w) view along direction y 1 = y 2 = y 2a = y 2b
F(s, t) =   R cos(s) R(sin(s) + 1) t   with s ∈ [0, 2π] and t ∈ -h 2 , h 2
where h is large enough to cover the entire desired length (h > h p cos(γ)).

Rotation between R 1 and R 2 will then be γ.

Translation between R 1 and R 2 will be given by h p /2 -h o .

3. Reference R 2a related to P 0 :

-origin P 0 -vector y 2a : equal to y 1 -vector x 2a : perpendicular to y 2a and contained in the plane defined by y 2a and C 0 (u) -vector z 2a : result of the vector product between x 2a and y 2a .

4. Reference R 2b related to P 1 : defined as R 2a for P 1 and C 1 (u)

Rotation between references R 2a and R 1 will be given by

θ 0 = arctan T 0 • z 1 T 0 • x 1 Similarly, between R 2b and R 1 we obtain θ 1 = arctan T 1 • z 1 T 1 • x 1
As translation between these references is h p /2, we can therefore calculate the passage matrices between R 1 and R 2a or R 2b .

We can thus express C 0 (u) and C 1 (u) in the references R 2a and R 2b : C 0 r 2a (u) and C 1 r 2b (u).

Resolution

The condition for tangency between the tool and the directrice C 0 (u) is translated by the unicity of the intersection between C 0 (u) and the ellipse Ell 0 (w) defined as the intersection of the tool and the plane (P 0 , x 2a , y 2a ). By defining tanEll 0 (w) as the vector tangent to the ellipse Ell 0 (w) we can express the unicity of the intersection by the fact that tanEll 0 (w) and C 0 (u) have the same slope in reference R 2a , at the point of intersection.

We can calculate the coefficients for sign µ, λ 0 and λ 1 such that:

If γ > 0, µ = 1, if γ < 0, µ = -1 where γ is the angle oriented (z 1 , z 2 ) If θ 0 > 0, λ 0 = 1, if θ 0 < 0, λ 0 = -1 where θ 0 is the angle oriented (z 1 , z 2a ) If θ 1 > 0, λ 1 = 1, if θ 1 < 0, λ 1 = -1 where θ 1 is the angle oriented (z 1 , z 2b )
Further, we can postulate

γ 0 = |γ| whence γ = µγ 0
The equation for Ell 0 (w) will be as follows, in the reference R 2a , (see Fig. 5):

Ell 0 (w) =   a 0 cos(w) + µL 0 R sin(w) + R 0   w ∈ [0, 2π] (1)
The vector tanEll 0 (w) tangent to Ell 0 (w) will be as follows, in the reference R 2a (see Fig. 5):

tanEll 0 (w) = ∂Ell 0 (w) ∂w ∂Ell 0 (w) ∂w =   -a 0 sin(w) R cos(w) 0   w ∈ [0, 2π] (2)
Ell 0 (w) The vector T 0 (u) tangent to C 0 (u), in the reference R 2a will be defined by

P 0 C 0 (u) tanEll 0 (w 0 ) M 0 x 2a y 2a R a 0 L 0 w 0 β 0 view along direction z 2a
T 0 (u) = ∂C 0 (u) ∂u ∂C 0 (u) ∂u u ∈ [0, 1]
and its slope β 0 will be given by

β 0 = arctan T 0 • y 2a T 0 • x 2a (3)
Further, we shall obtain (see Fig. 4):

tan(γ 0 ) = L 0 cos(λ 0 θ 0 ) h 0 -µλ 0 L 0 sin(λ 0 θ 0 ) whence:

L 0 = h 0 tan(γ 0 ) cos(λ 0 θ 0 ) + µλ 0 tan(γ 0 ) sin(λ 0 θ 0 ) (4)
and (see Fig. 4):

cos(λ 0 θ 0 -µλ 0 γ 0 ) = R a 0 whence: a 0 = R cos(λ 0 (θ 0 -µγ 0 )) (5)
Similarly, in the reference R 2b we have: The equation for Ell 1 (w):

Ell 1 (w) =   a 1 cos(w) -µL 1 R sin(w) + R 0   w ∈ [0, 2π] (6)
The vector tanEll 1 (w) tangent to Ell 1 (w):

tanEll 1 (w) = ∂Ell 1 (w) ∂(w) ∂Ell 1 (w) ∂(w) =   -a 1 sin(w) R cos(w) 0   w ∈ [0, 2π] (7)
The vector T 1 (u) tangent to C 1 (u):

T 1 (u) = ∂C 1 (u) ∂u ∂C 1 (u) ∂u u ∈ [0, 1]
whence its slope

β 1 = arctan T 1 • y 2b T 1 • x 2b (8)
In the same way as in R 2a reference, we obtain:

L 1 = h 1 tan(γ 0 ) cos(λ 1 θ 1 ) + µλ 1 tan(γ 0 ) sin(λ 1 θ 1 ) (9)
and:

a 1 = R cos(λ 1 (θ 1 -µγ 0 )) (10)
Using ( 4) and ( 5) in ( 1) and (2), we obtain the equations at the point of contact M 0 :

The intersection at M 0 will be translated by

C 0 r 2a (u 0 ) • x 2a = Ell 0 (w 0 ) • x 2a and C 0 r 2a (u 0 ) • y 2a = Ell 0 (w 0 ) • y 2a
That is,

C 0 r 2a (u 0 ) • x 2a = R cos(λ 0 (θ 0 -µγ 0 )) cos(w 0 ) + µ h 0 tan(γ 0 ) cos(λ 0 θ 0 ) + µλ 0 tan(γ 0 ) sin(λ 0 θ 0 ) ( 11
)
and

C 0 r 2a (u 0 ) • y 2a = R sin(w 0 ) + R (12)
The equality of slopes between tanEll 0 (w) and C 0 r 2a (u) will be translated by

β 0 = tanEll 0 • y 2a tanEll 0 • x 2a
That is

β 0 = -cos(λ 0 (θ 0 -µγ 0 )) tan(w 0 ) (13)
Similarly, we shall obtain the following in reference R 2b :

C 1 r 2b (u 1 ) • x 2b = R cos(λ 1 (θ 1 -µγ 0 )) cos(w 1 ) -µ h 1 tan(γ 0 ) cos(λ 1 θ 1 ) + µλ 1 tan(γ 0 ) sin(λ 1 θ 1 ) (14) C 1 r 2b (u 1 ) • y 2b = R sin(w 1 ) + R ( 15
)
β 1 = -cos(λ 1 (θ 1 -µγ 1 )) tan(w 1 ) (16)
Further, we know that

h 0 + h 1 = h p (17)
Whence a system of 7 equations with 7 unknown values (u 0 , u 1 , w 0 , w 1 , h 0 , h 1 and γ) that can be resolved numerically in a sufficiently short time.

3 Error Calculation and Analysis

Error Calculation Methodology

The errors on both sides of point M 2 are calculated separately.

Comment: the iteration increment between two successive rules will then be considered small enough for the two successive positions of the tool to be included within the envelope of tool positions. Otherwise, the error would need to be calculated in relation to this envelope, which would require an excessively long calculation time.

The calculation is carried out by postulating

S(u, v) = F(s i , t i ) + ε n T (s i , t i )
where F(s i , t i ) is the point of the tool where the calculation is carried out, s i and t i are values for parameters at that point, n T (s i , t i ) the unit vector perpendicular to the tool at that point, and ε the value sought for. This equation can thus be reduced to a system of three equations and three unknown values: u, v and ε. In this equation, s i and t i are taken between the limits of the parametric interval for calculation determined by the contact points. We can then choose the values for parameters at contact points M 0 , M 1 , and M 2 as limits for our calculation interval.

The resolution algorithm adopted will be a dichotomial algorithm with variation in the two parametric directions s and t.

Between M 2 and M 0 , we take the following as our calculation interval:

s 2 < s i < s 0 t 2 < t i < t 0
While between M 2 and M 1 , we take the following as our calculation interval:

s 1 < s i < s 2 t 1 < t i < t 2

Analysis Principle

Given that we cannot carry out an analytical error calculation, the error calculation study is made through a range of representative examples. For this purpose, we take into account the fact that the main parameters having an influence on error are:

-The tool radius R.

-The angle α made by the two perpendiculars N 0 and N 1 . -The length of the rule h p . We will study the specific influence of each parameter on the error made through tests carried out on a hyperbolic paraboloid, allowing us to change each parameter independently.

The hyperbolic paraboloid we studied was based on the two following directrices:

C 0 (u) =   k u k tan(α 0 ) u -1 2 hp 2   and C 1 (u) =   k u -k tan(α 0 ) u -1 2 -hp 2  
whence the parametered surface

S(u, v) =   k u -1 2 k tan(α 0 )(4uv -2u -2v + 1) 1 2 h p -v h p  
where k is a parameter for scale, we shall take k = h p .

In order to validate this approach, we compared the error curves obtained with any ruled surface (that defined in (18), see Fig. 6) and a hyperbolic paraboloid (see Fig. 7). We can observe qualitatively similar curve shapes. The setting we suggest for the tool would appear to allow for a considerable reduction in interference between the tool and the surface.

Further, the calculation time needed to determine setting parameters and error is no greater than a few seconds, which leads us to think that implementation in CAD/CAM type software is feasible.

As an example, we carried out comparative analysis for any ruled surface:

x 1 z 1 x 1 z 1 Figure 8: Error localisation S(u, v) =   30u 2 + 10u + 30 + 20vu 2 -10vu -10v 100 √ u -50 -100v √ u + 90v -60v sin 1 2 πu 70 -30u 2 -130v + 50vu 2   (18)
For this surface, we chose u p = 0.4 and R = 40. With standard setting with error distribution, max interference will be ε = 4.880mm.

With the optimum setting we suggest max interference will be ε = 0.181mm.

Figure 8 allows error localisation in each case to be appreciated qualitatively.

Figure 8 also allows another advantage of the suggested model to be highlighted, which is to respect the directrices of the ruled surface which can prove to be extremely useful in a large number of cases, particularly within the perspective of a circumferential machining method for any form of surface using piecewise approximation as proposed in [START_REF] Elber | 5 Axis free-form surface milling using piecewise ruled surface approximation[END_REF]. We can thus envisage a side milling methodology for complex surfaces, which is particularly valid for convex and double curved surfaces, allowing for a considerable gain in time in relation to conventional end machining techniques.

Further, this model enabled us to develop a computation algorithm for the optimum tool radius for a given case. The major difficulty encountered was in optimising the algorithm to obtain this radius with a minimum number of iterations. We chose the ratio of error for standard setting over error for optimised setting as the indicator for performance of our setting as compared with standard positioning. All the results obtained hitherto show that the model we suggest allow errors to be reduced by a ratio of at least 10 to !.

Further, for the lowest performance indicator obtained (e o = 45), it can be observed that the gain is greater from a quantitative point of view.

Figure 4 :

 4 Figure 4: Definition of references.

Figure 5 :

 5 Figure 5: Ellipse Ell 0 (w)

ErrorFigure 6 :ErrorFigure 7 :

 67 Figure 6: Error curves for any ruled surface

Figure 9 :

 9 Figure 9: Performance ratings in relation to h p .

  Figures 9 and 10 show the variation of this indicator in relation to R for different values of Oto on the one hand and hp o n the other hand.

Using the model we suggest, for the same radius R = 10 mm, we obtained ε = 0.22 mm, i.e. a performance ratio of (0.585 × 2)/0.22 = 5.3.


Sonia Cafieri

Pierre Hansen

Nenad Mladenović

Variable Neighborhood Search for edge-ratio network clustering

Edge-ratio clustering was introduced in [Cafieri et al., Phys.Rev. E 81(2):026105, 2010], as a criterion for optimal graph bipartitioning in hierarchical divisive algorithms for cluster identification in networks. Exact algorithms to perform bipartitioning maximizing the edge-ratio were shown to be too time consuming to be applied to large datasets. In this paper, we present a Variable Neighborood Search (VNS)-based heuristic for hierarchical divisive edge ratio network clustering. We give a full description including the structure of some algorithmic procedures which are used to implement the main steps of the heuristic. Computational results show that the proposed algorithm is very efficient in terms of quality of the bipartitions, moreover the computing time is much smaller than that one for exact algorithms.

Introduction

Network analysis is successfully used in the study of complex systems in a variety of domains, where a network representation and the study of its topological features help to better understand some characteristics of the system under consideration. Prominent examples include social networks, describing individuals and their interactions and relationships, telecommunication networks, such as the World Wide Web, transportation networks, biological networks, and many more. A detailed introduction to networks has recently been given by Newman [START_REF] Newman | Networks: an introduction[END_REF]. A mathematical representation of a network is obtained using a graph G = (V, E), where vertices (in the set V ) are associated to the entities under study and the edges (in the set E) joining pairs of vertices correspond to relationships among the entities.

Network clustering represents a topic of particular interest in network analysis. It consists in detecting subsets of vertices which are more densely linked compared to the rest of the graph, called modules or clusters or communities. This is of much interest in the study of complex systems as many of such systems are characterized by a modular structure. For example, a community in a social network can be constituted by individuals sharing a common interest or location, in a biological network by entities with a common function. The reader is referred to Fortunato [START_REF] Fortunato | Community detection in graphs[END_REF] for a recent survey of this research domain.

Network clustering is generally based on a definition of community or on a criterion to evaluate a partition found. A clustering criterion can be used as an objective function in a combinatorial optimization problem, whose solution gives an optimal partition for the considered network. Several alternatives have been proposed. Among clustering criteria, the most used is modularity [START_REF] Newman | Finding and evaluating community structure in networks[END_REF], based on the idea of comparing the fraction of edges falling within communities to the expected fraction of such edges. Two of the most known definitions of community were proposed in [START_REF] Radicchi | Defining and identifying comminities in networks[END_REF] and give conditions defining communities in a strong and a weak sense, respectively. Building on the most significant condition, the weak condition, in [START_REF] Cafieri | Edge ratio and community structure in networks[END_REF] a new criterion was introduced for a bipartition to be optimal, called the edge-ratio criterion. This criterion, as well as the weak condition, are recalled below. The edge-ratio criterion is maximized to obtain successive bipartitions in a hierarchical divisive algorithm. In [START_REF] Cafieri | Edge ratio and community structure in networks[END_REF] the optimal bipartitioning problem was solved exactly at each iteration of the divisive algorithm. It turns out that this exact algorithm is very time consuming even for small or medium size datasets (i.e., with a few hundred of vertices). Therefore, heuristic approaches appear to be the most adequate to solve in reasonable time the bipartitioning problems. In this paper, we propose a Variable Neighborhood Search heuristic to solve the bipartitioning problems. Such heuristic was presented in [START_REF] Cafieri | Edge-ratio network clustering by variable neighborhood search[END_REF] somewhat informally. In the present paper we give a full mathematical description of this heuristic, including the structure of some algorithmic procedures which are used to implement the main steps of the heuristic.

Let us consider a graph G = (V, E) with V set of vertices, with cardinality n, and E set of edges, with cardinality m. Let S ⊆ V be a subset of vertices. Then the degree k i of a vertex i belonging to V (i.e., the number of its neighbors) can be separated into two components k in i (S) and k out i (S), which represent respectively the number of neighbors of i inside S and the number of neighbors of i outside S. Radicchi et al. [START_REF] Radicchi | Defining and identifying comminities in networks[END_REF] give the following definition for a set of vertices S forming a community in the weak sense:

i∈S k in i (S) > i∈S k out i (S)
That is, S is a community in the weak sense if and only if the sum of all degrees within S is larger than the sum of all degrees joining S to the rest of the graph. As this condition can be used ( [START_REF] Radicchi | Defining and identifying comminities in networks[END_REF]) as a local stopping criterion in hierarchical clustering, it has been used by Wang et al. [START_REF] Wang | Remarks on network community properties[END_REF] to define a community S indivisible if there is no bipartition, (S 1 , S 2 ) of S, such that both S 1 and S 2 satisfy the weak condition.

In [START_REF] Cafieri | Edge ratio and community structure in networks[END_REF] the definition of community in the weak sense is extended into a criterion for a bipartition to be optimal: one seeks to maximize the minimum, for both classes S 1 and S 2 of the bipartition of S (such that

S 1 ∪ S 2 = S, S 1 ∩ S 2 = ∅, S 1 , S 2 = ∅)
, of the ratio of inner edges to cut edges:

f (S 1 , S 2 ) = max S 1 ,S 2 ⊂V min (r(S 1 ), r(S 2 )) ,
where for a subgraph S the ratio r(S) is given by:

r(S) = i∈S k in i (S)/ i∈S k out i (S).
This definition is strenghten to be used in a hierarchical divisive clustering algorithm, by introducing a parameter α and by quantifying how much the number of inner edges is larger than the number of cut edges. Thus, the condition becomes:

i∈S k in i (S) ≥ α i∈S k out i (S). (1)
In case of equality, the coefficient α is equal to the ratio of twice the number of edges within the community S divided by the number of edges within the cut of that community. It is called the edge ratio [START_REF] Cafieri | Edge ratio and community structure in networks[END_REF]. Then, hierarchical divisive clustering algorithm, where at each step a bipartitioning of the current subgraph is done, one can then seek the maximum value of α for which the subgraph to be bipartitioned will be divisible. So, the problem of detecting indivisible communities is also solved. The crucial step of bipartitioning in a hierarchical divisive clustering algorithm is also the most computationally demanding. The VNS based heuristic presented in this paper is aimed at performing efficiently this algorithmic step.

The paper is organized as follows. In Sect. 2 we recall the structure of a hierarchical divisive clustering algorithm and introduce VNS to perform its main computational step. Sect. 3 presents the details of the algorithmic structure of some subroutines which are used to implement the main steps of the heuristic. A few results about algorithmic complexity are also discussed. In Sect. 4 the results of computational experiments performed on a set of instances from the literature are presented and discussed. Sect. 5 concludes the paper.

2 Divisive hierarchical clustering using VNS

Basic notation

Let us consider an unweighted graph G(V, E) with vertex set V = {v 1 , . . . , v n } and edge set E = {e 1 , . . . , e m }. Let k i denotes the degree of vertex v i ∈ V . Let A = a ij be the corresponding adjacency matrix (a ij ∈ {0, 1}) and H = h ij be the adjacency list, i.e., raw i contains the list h ij , j = 1, . . . , d i of vertices that are adjacent to vertex v i .

Divisive hierarchical clustering

Hierarchical divisive clustering heuristics (see, e.g. [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF]) proceed, in a top-to-bottom approach, from an initial partition of the graph containing all its n vertices, iteratively bipartitioning a subgraph, until a partition into n clusters is obtained, or the value of a chosen objective function is not improved anymore. At each step of a divisive clustering algorithm, one must select a cluster among the ones obtained at the previous hierarchy levels and divide it into two. At the first step, the first bipartition is done on V . Let ℓ in {1, . . . , n} be the index of the subgraph S ℓ selected to be bipartitioned during a step of the divisive algorithm. After a bipartition into subsets of vertices S 1 and S 2 , we obtain S ℓ = S 1 and S ℓ+1 = S 2 . Then ℓ is increased by one and all process repeated. In the present study, the objective function is based on the edge-ratio criterion [START_REF] Cafieri | Edge ratio and community structure in networks[END_REF]. So, for each bipartition, the objective function value f (S 1 , S 2 ) is computed. If it is greater than or equal to one, i.e., the subgraph that one seeks to bipartition is divisible, then the obtained subgraphs are added to the list P of subgraphs to be considered for bipartition. The corresponding pseudo-code is given in Algorithm 1. Notice that bipartitioning is done by Variable Neighbourhood Search, so an appropriate function BVNS is called to perform the bipartitioning steps. It is described in details in Sect. 3.

Function DHC-VNS (m, n, A, P )

1 ℓ = 1, S ℓ = {v 1 , . . . , v n }, P = ∅ 2 while |S ℓ | > 1 do 3 S 1 , S 2 ← BVNS (S ℓ ) // Find bipartition of S ℓ according to f (S 1 , S 2 ); 4 if f ≥ 1 then 5 P = P ∪ arg min f (S 1 , S 2 ) 6 S ℓ ← S 1 ; S ℓ+1 ← S 2 7 ℓ ← ℓ + 1 8
Among {S 1 , . . . , S ℓ }, select cluster with the largest cardinality 9 Exchange its position with S ℓ Algorithm 1: Divisive hierarchical clustering algorithm with VNS for bipartitioning

Variable Neighborhood Search for edge-ratio clustering

Variable Neighbourhood Search (VNS) is a metaheuristic aimed at solving combinatorial and global optimization problems. It is based on the idea of a systematic change of neighbourhood combined with a local search to escape the current local optimum. The main ingredients of such heuristic are thus given by a local search procedure, shaking, i.e., a procedure to perturb the current solution, and a procedure to update the current solution. The reader is referred to [START_REF] Mladenović | Variable neighbourhood search[END_REF][START_REF] Hansen | Variable neighbourhood search: Principles and applications[END_REF][START_REF] Hansen | Variable neighbourhood search[END_REF][START_REF] Hansen | Variable neighbourhood search[END_REF][START_REF] Hansen | Variable neighbourhood search: methods and applications[END_REF] for an introduction of VNS and its main applications.

In this section, a detailed description of the proposed Variable Neighborhood Search is given, together with pseudo-codes for its main procedures.

Initialization

The VNS for bipartitioning starts with a random partition of the current vertex set S = S 1 ∪S 2 , S 1 ∩ S 2 = ∅, where the cardinality of the first set is n

2 (|S 1 | = [ n 2 ]
). In that way, two subgraphs

G 1 = (S 1 , E 1 ) and G 2 = (S 2 , E 2 ) are obtained. A third set of vertices S 3 = V \ {S 1 ∪ S 2 }
is also handled through the execution of the algorithm. Obviously, at the first iteration S 3 = ∅. Let P ℓ denotes all possible partitions of vertex set S ℓ , where ℓ represents the hierarchy level of the divisive hierarchical clustering algorithm. In order to simplify notation, in the following we will omit the index ℓ.

Data structure and objective function

Let S be a set of vertices to be bipartitioned into S 1 and S 2 . A bipartition is represented in the computer memory in the following way:

x = {x 1 , . . . , x |S| } -indices of vertices that belong to S δ x j = 1 if x j ∈ S 1 0 if x j ∈ S 2
and the following variables are used to update the objective function values:

n 1 -
k rest i = k i -k i (S), v i ∈ S -"rest" degree of vertex v i ,
where k i (S) denotes the degree of the vertex in the subgraph of G induced from vertices of S;

f 1 = 2 • m 1 c 1 -ratio for subgraph G 1 ; f 2 = 2 • m 2 c 2 -ratio for subgraph G 2 ; f = min{f 1 , f 2 } -objective function value.
Note that all values from above can be computed if a partition (S 1 , S 2 , S 3 ) is known (S j ∈ P). Therefore, instead of n 1 we could in fact put n 1 (S 1 ), instead of m 2 we could put m 2 (S 2 ), etc. Also note that the degree of each vertex v can be presented as

k v = k in v + k cut v + k rest v
with respect to the current partition of V in three subsets.

At each step of the divisive hierarchical clustering algorithm, the problem is to find a bipartition of the current S into S 1 and S 2 such that the largest value between the objective functions associated to the corresponding subgraphs, f 1 and f 2 , is as small as possible. The objective function of the edge-ratio clustering at each iteration of the divisive algorithm can be presented as follows

f (S 1 , S 2 ) = max S 1 ,S 2 ∈P min{f 1 (S 1 ), f 2 (S 2 )} (2)
where S 1 ∪ S 2 = S and S 1 ∩ S 2 = ∅.

Proposition 1 Given a graph G(V, E) and the bipartition of its vertex set V into sets S 1 and S 2 , computing the objective function value (2) is in O(|E|).

Proof. In order to compute f , we need to find m 1 , m It consists in taking some entity v from S = S 1 ∪ S 2 and changing its allocation: if v ∈ S 1 then, after the allocation move, it will belong to S 2 and vice-versa.

The vertex whose allocation change improves the most the objective function ( 2) is selected and all variables updated accordingly. The datailed algorithms are given below.

Move.

In Algorithm Move(v, f, δ) input values are the index of entity v that will change its membership, the current objective function value f and the current solution δ. We also assume that values of m 1 , m 2 , c 1 and c 2 are known as well.

Function Move(v, f, δ)

1 if δ v then 2 mm 1 ← m 1 -k in v ; mm 2 ← m 2 + k cut v 3 cc 1 ← c 1 + k in v -k cut v -k rest v 4 cc 2 ← c 2 + k in v -k cut v + k rest v else 5 mm 1 ← m 1 + k cut v ; mm 2 ← m 2 -k in v 6 cc 1 ← c 1 + k in v -k cut v + k rest v 7 cc 2 ← c 2 + k in v -k cut v -k rest v 8 f 1 ← 2•mm 1 cc 1 ; f 2 ← 2•mm 2 cc 2 ; f ← min{f 1 , f 2 }
Algorithm 2: Allocation move Notice that mm 1 , mm 2 , cc 1 and cc 2 are used to temporally denote the values of m 1 , m 2 , c 1 and c 2 ; these values will be actually changed once the move to be performed is selected.

Proposition 2 Let v be a vertex that moves from one subset to another. If v ∈ S 1 then the new objective function value is

f new = min{ 2 • (m 1 -k in v ) c 1 + k in v -k cut v -k rest v , 2 • (m 2 + k cut v ) c 2 + k in v -k cut v + k rest v }. If v ∈ S 2 then f new = min{ 2 • (m 1 + k cut v ) c 1 + k in v -k cut v + k rest v , 2 • (m 2 -k in v ) c 2 + k in v -k cut v -k rest v }.
Proof. A new value of the objective function, by definition is

f new = min{f new 1 , f new 2 } = min{ 2m new 1 c new 1 , 2m new 2 c new 2 }. (3)
If v ∈ S 1 changes its membership to S 2 , then we have:

m new 1 = m 1 -k in v ( 4
)
m new 2 = m 2 + k cut v (5)
because the new number of inner edges m 1 in S 1 is reduced by k in v (the inner degree of v) (see ( 4)), and the number of edges in S 2 is augmented by the cut degree of v (see [START_REF] Girvan | Community structure in social and biological networks[END_REF]). The numbers of cut edges c 1 and c 2 of S 1 and S 2 are updated as follows:

c new 1 = c 1 + k in v -k cut v -k rest v (6) c new 2 = c 2 + k in v -k cut v + k rest v (7)
Indeed, since v moves to S 2 , the number of cut edges c new 1 of S 1 is increased by its inner degree in S 1 (+k in v ); also, it is reduced by k cut v , since its cut degree contributes to inner degree of S 2 ; finally, all vertices from S \ (S 1 ∪ S 2 ) connected with v after the move to S 2 do not produce cut edges of S 1 anymore (-k rest v ). Observe that c new 1 = c new 2 only if S = V , i.e., at the first hierarchy level ℓ = 1, where subset S 3 = ∅ and therefore the rest degrees are equal to 0. Similarly for [START_REF] Hansen | Variable neighbourhood search: Principles and applications[END_REF].

If vertex v belongs to S 2 and moves to S 1 , then new numbers of inner and cut edges in S 1 and S 2 are

m new 1 = m 1 + k cut v (8) m new 2 = m 2 -k in v (9) c new 1 = c 1 + k in v -k cut v + k rest v ( 10
)
c new 2 = c 2 + k in v -k cut v -k rest v (11)
Substituting ( 6) -( 11) into equations (3) we get the result. ⋄

Update.

Once the entity v improving the most the objective function value is selected, a move is done and the current solution as well as all the variables that we keep during the search are updated.

In order to get a new value improving the objective function, both f new 1 and f new 2 have to be larger than f , since the minimum of these two values should be larger than f . Let us first assume that v ∈ S 1 is moved to S 2 . Besides updating values with formulas ( 4) -( 11), we Function Update(v, f, δ)

1 if δ v then 2 n 1 ← n 1 -1; n 2 ← n 2 + 1 3 m 1 ← m 1 -k in v ; m 2 ← m 2 + k cut v 4 c 1 ← c 1 + k in v -k cut v -k rest v 5 c 2 ← c 2 + k in v -k cut v + k rest v else 6 n 2 ← n 2 -1; n 1 ← n 1 + 1 7 m 1 ← m 1 + k cut v ; m 2 ← m 2 -k in v 8 c 1 ← c 1 + k in v -k cut v + k rest v 9 c 2 ← c 2 + k in v -k cut v -k rest v 10 y ← k cut v ; k cut v ← k in v ; k in v ← y 11 for j = 1, k v do 12 r ← h ij 13 if (δ r ∧ δ v ) ∨ (¬δ r ∧ ¬δ v ) then 14 k cut r ← k cut r + 1; k in r ← k in r -1 else 15 k cut r ← k cut r -1; k in r ← k in r + 1 16 f 1 ← 2•m 1 c 1 ; f 2 ← 2•m 2 c 2 ; f ← min{f 1 , f 2 } 17 δ v ← ¬δ v ;

Algorithm 3: Updating values after allocation move

In line 1 we check if v ∈ S 1 . If so, all variables are updated: in line 2 the cardinality of S 1 is reduced by one (n 1 ← n 1 -1), and cardinality of S 2 increased by one (n 2 ← n 2 + 1). The total number of inner and cut degrees of both sets (S 1 and S 2 ) are updated in lines 3, 4 and 5 respectively. From step 11 to 15 the degrees of vertices adjacent to v are updated. Similar steps are repeated if v ∈ S 2 . Finally, the new objective function is computed, and membership of vertex v flipped.

The following proposition can be easily deduced looking at the given pseudo-code.

Proposition 4

The worst case complexity of algorithm Update is in O(k max ), where k max is the maximum degree of G.

Note that, comparing the result of Proposition 2 with the result of Proposition 1, we see that the number of operations is decreased by order of mk max by updating.

Local search algorithm.

A local search, or steepest descent procedure, is built on the idea of moving vertices from one subset to another until an improvement in the objective function value, i.e., until f new > f , is obtained. The vertex whose change of membership produces the maximum improvement is chosen to make a new partition. Such partition corresponds to a local maximum of the optimization problem.

Algorithm Best-Impr-LS(x, f opt , δ) describes the steps of the proposed Allocation local search.

Function Best-Impr-LS(x, f opt , δ)

1 f opt ← ∞; improve ← true 2 while Improve do 3 Improve ← false 4 for i = 1, n do 5 v ← x i 6 if k cut v > 0 then 7 Move(v, f, δ) 8 if f < f opt then 9 f opt ← f ; w ← v 10 improve ← true 11 if ¬improve return 12 Update(w, f, δ) Algorithm 4: Best improvement Allocation (re-assignment) local search
The following property holds.

Proposition 5

The time complexity of one iteration of the local search algorithm is O(n).

Proof. The proposed local search is based on the move procedure, which is repeated n times. The results follows from the observation that the procedure move is O(1). ⋄

Shaking

Since the solution of bipartitioning S is represented with a binary vector δ(x j ), the natural way to introduce a distance into the solution space is to use the Hamming distance. That is, we say that the two solutions x and y are at distance k if their Hamming distance is equal to k, i.e., the minimum number of changes of zeros to ones (or reciprocally) required to change x into y is equal to k. We build the shaking procedure on this concept, perturbing the current solution x to a new vector which has Hamming distance k from it. This means that k zeros or ones are switched to their opposite value. In other words, k entities change their original cluster. Therefore, our Shake(k, δ) procedure works as described in Algorithm 5.

Function Shake(k, δ)

1 ℓ ← 1 2 while ℓ < k do 3 r ← 1 + (|S| -ℓ) • Rand 4 v ← x r 5
Update(v, f, δ)

6 ℓ ← ℓ + 1;
Algorithm 5: Steps of Shaking operator

Basic VNS for edge-ratio clustering

The basic VNS heuristic for edge-ratio clustering is straightforward, once the local search and the shaking procedures have been introduced. It is summarized in Algorithm 6. Note that BVNS is called at each iteration of the divisive hierarchical clustering Algorithm 1 to perform bipartitioning.

Function BVNS (δ, k max , t max )

1 Initialization(δ) 2 repeat 3 k ← 1 4 repeat 5 δ ′ ← Shake(k, δ) /*

Computational results and analysis

The proposed VNS heuristic for edge-ratio clustering is validated on a set of 11 problems from the literature, corresponding to various real-world applications. They are: Zachary's karate club dataset [START_REF] Zachary | An information flow model for conflict and fission in small group[END_REF], describing friendship relationships in a karate club; Lusseau's dolphins dataset [START_REF] Lusseau | The bottlenose dolphin community of doubtful sound features a large proportion of longlasting associations. can geographic isolation explain this unique trait?[END_REF], describing communications between dolphins; Hugo's Les Misérables dataset [START_REF] Knuth | The Stanford GraphBase: A Platform for Combinatorial Computing[END_REF], describing characters in Victor Hugo's masterpiece and their interactions; Krebs' political books dataset [START_REF] Krebs | [END_REF]; a dataset representing the schedule of football games between American college teams [START_REF] Girvan | Community structure in social and biological networks[END_REF]; a dataset dealing with connections between US airports [19]; a dataset on a coauthorship network of scientists working on network theory and experiment [START_REF] Newman | Finding community structure in networks using the eigenvectors of matrices[END_REF]; a network describing electronic circuits [START_REF] Alon | [END_REF]; a network representing e-mail interchanges between members of a university [START_REF] Guimerà | Self-similar community structure in a network of human interactions[END_REF]; a network giving the topology of the Western States Power Grid of the United States [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF] and a network of authors collaborations [19]. All considered networks are undirected and unweighted, without loops.

We compare the results obtained using the VNS proposed in this paper for edge-ratio bipartitioning in a divisive hierarchical algorithm (Algorithm 1) with the results obtained solving, at each iteration of the divisive algorithm, the bipartition problem by an exact algorithm. An exact solution of the bipartitioning problem based on the edge-ratio criterion was proposed in [START_REF] Cafieri | Edge ratio and community structure in networks[END_REF]. In that paper, a mathematical programming formulation is proposed where the edge-ratio α is maximized subject to non-linear and non-convex constraints, coming from the products of α with binary variables. An exact algorithm is then obtained fixing the value of α to solve a linear program in 0-1 variables, and then proceeding with a dichotomous search on the values of α. Specifically, one strarts with a value of α equal to 1. If there is no feasible solution for that value, the network is indivisible and one stops, otherwise the value of α is doubled and feasibility checked until a value is attained for which the weak condition cannot be satisfied. This gives an upper bound ᾱ and the previous value of α gives a lower bound α. Then the dichotomous search proceeds by considering the mid value of the interval [α, ᾱ] (see [START_REF] Cafieri | Edge ratio and community structure in networks[END_REF]). It turns out that this exact algorithm is very time consuming even for small or medium size datasets (i.e., with a few hundred of vertices).

We implemented BVNS with the value of k max set to min{50, n/2}. Tests were run on AMD Opteron 2 GHz CPU, 128 GB RAM. Results are shown in Table 1. Networks are listed together with their number of vertices n and number of edges m. We report the cpu time (seconds) to obtain clustering solutions by a divisive hierarchical algorithm using the VNS and the above exact algorithm for the bipartitioning problem solved at each hierarchical level. problem at each hierarchical level significantly faster than exact bipartitioning, reducing the computational time up to 7 orders of magnitude. This allows us to solve larger problems, increasing the number of problems that we can solve using the exact method. In terms of quality of the solution, the proposed VNS also appears to be efficient. In the most of the cases, we obtain the same partition using VNS and the exact bipartitioning method. We give in Figures 1 and2 two examples of obtained partitions, and the corresponding dendrogram illustrating the hierarchical divisive algorithm.

Figure 1 :Figure 2 :

 12 Figure 1: Partition and dendrogram obtained for dataset polbooks

  number of vertices in set S 1 ; n 2 -number of vertices in set S 2 ; n 3 -number of vertices in set S 3 ; m 1 -number of edges in subgraph G 1 ; m 2 -number of edges in subgraph G 2 ; c 1 -number of cut edges with one end point in S 1 and another in V \ S 1 = S 2 ∪ S 3 ; c 2 -number of cut edges with one end point in S 2 and another in V \ S 2 = S 1 ∪ S 3 ;

	k in

i -inner degree of vertex v i (both end points belong to either S 1 or S 2 ); k cut i -cut (outer) degree of vertex v i (number of links of v i with vertices from different sets);

  [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF] , c 1 and c 2 . For each vertex v i we check if h ij belongs to S 1 or S 2 (i.e., we check if δ(h ij ) = true or false). If both δ i and δ(h i,j ) have the same value, the corresponding inner degree (m 1 or m 2 ) is increased by one. Otherwise, cut degree (c 1 or c 2 ) is increased by one. Since the number of elements in H is m = |E| we need O(m) calculations. ⋄3.3 Allocation local searchLocal search in the proposed VNS is based on the selection of the entity (vertex) whose associated local change (move) gives the best improvement of the objective function value, and on the update of all entities involved in the problem after that move. Assume that any solution (S 1 , S 2 , S 3 ) is known. We use an allocation move as a local change.

  Basic VNS for edge-ratio clusteringIn Algorithm 6, k max and t max represent respectively the maximum number of explored neighborhoods and the maximum cpu time used in the stopping criterion.

			Shaking */
	6	δ ′′ ← Best-Impr-LS(δ ′ , f )	/* Local search */
	7	k ← k + 1	/* Next neighborhood */
	8	if f (δ ′′ ) < f (δ) then
	9	δ ← δ ′′ ; k ← 1	/* Make a move */
		until k = k max
	10	t ← CpuTime()
		until t > t max
		Algorithm 6:

Table 1 :

 1 Comparison of results obtained using VNS and an exact algorithm for bipartition according to the edge-ratio criterion, on datasets from the literature It appears that the proposed VNS for edge-ratio is very efficient in terms of quality of ob-

	dataset	n	m	Time VNS (sec.) Time exact (sec.)
	karate	34	78	0.0	62.10
	dolphins	62	159	9.99e-3	172.2
	les miserables	77	254	1.99e-2	283.49
	political books	105	441	1.99e-2	716.45
	football	115	613	3.99e-2	11780.79
	Usair97	332	2126	0.62	3752906.94
	netscience main	379	914	1.14	2347.24
	s838	512	819	0.64	12612.53
	email	1133 5452	6.84	-
	power	4941 6594	29.80	-
	erdos02	6927 11850	35.50	-

ConclusionsWe proposed a Variable Neighborhood Search heuristic for solving the bipartitioning problem in a hierarchical divisive graph clustering algorithm, according to the recently introduced edge-ratio criterion. Neighborhoods to be used in the VNS are defined using the Hamming distance. We describe the proposed local search and shaking procedures, the latter used to perform systematic perturbations of the current incumbent solution. We give a full description including the pseudo-codes of the procedures used to implement the main steps of the heuristic.

also need to update the number of vertices of S 1 (reduced by 1), and the cardinality of S 2 (increased by 1):

In addition, inner and cut degrees of vertices that are adjacent to the moved vertex v have to be updated.

Proposition 3 Let v be a vertex that changes its membership and let u be any adjacent vertex to v ((v, u) ∈ E). Then values of new inner and cut degrees of u are:

Proof. As mentioned earlier, u j = h vj , j = 1, . . . , k v indicate all adjacent vertices of v. There are k v of them since k v is the degree of v. Let us denote with kin u and kcut u a new values of inner and cut degrees (after the vertex v is moved from one to another subset of vertices), respectively. Then four different cases can occur. Case 1. v ∈ S 1 and u j ∈ S 1 . Then kin u j = k in u j -1 and kcut

Case 2. v ∈ S 1 and u j ∈ S 2 . Then kin u j = k in u j + 1 and kcut

Case 4. v ∈ S 2 and u j ∈ S 2 . Then kin u j = k in u j -1 and kcut

So, cases 1 and 4 give the same outcome in updating degrees of adjacent vertices to v, and the same holds for cases 2 and 3. Thus, we get the result. ⋄ Algorithm Update(i, f, δ) summarizes the above updating procedure. Input values are given by vertex v (that changes its membership), a solution δ and the objective function value f . Notice that all other values mentioned before (n 1 , n 2 , m 1 , etc.) are also input values, however we do not put them in the input list for the sake of readability.

Computational results show that VNS, compared to exact bipartitioning, gives good quality results significantly reducing the computational time.