Tawsif Gokhool

Renato Martins

Patrick Rives

Noëla Despré

A Compact Spherical RGBD Keyframe-based Representation

This paper proposes an environmental representation approach based on hybrid metric and topological maps as a key component for mobile robot navigation. Focus is made on an ego-centric pose graph structure by the use of Keyframes to capture the local properties of the scene. With the aim of reducing data redundancy, suppress sensor noise whilst maintaining a dense compact representation of the environment, neighbouring augmented spheres are fused in a single representation. To this end, an uncertainty error model propagation is formulated for outlier rejection and data fusion, enhanced with the notion of landmark stability over time. Finally, our algorithm is tested thoroughly on a newly developed wide angle 360 0 field of view (FOV) spherical sensor where improvements such as trajectory drift, compactness and reduced tracking error are demonstrated.

I. INTRODUCTION

Visual mapping is a required capability for autonomous robots and a key component for long term navigation and localisation. Inherent limitations of GPS in terms of dropping accuracy in densely populated urban environments or lack of observability in indoor settings have propelled the applications of visual mapping techniques. Moreover, the rich content provided by vision sensors maximises the description of the environment.

A metric map emerges from a locally defined reference frame and all information acquired along the trajectory is then represented in relation to this reference frame. Hence, localisation in a metric map is represented by a pose obtained from frame to frame odometry. The benefit of this approach is its precision at the local level allowing precise navigation. On the other hand, this representation is prone to drift phenomena which becomes significant over extended trajectories. Furthermore, the amount of information that is accumulated introduces considerable problems in memory management due to limited capacity.

By contrast, a topological map provides a good trade-off solution to the above-mentioned issues. Distinct places of the environment are defined by nodes and joined together by edges representing the accessibility between places generating a graph. This representation not only brings about a good abstraction level of the environment but common tasks such as homing, navigation, exploration and path planning become more efficient. However, the lack of metric information may render some of these tasks precision deficient as only rough *This work is funded by EADS/Airbus under contract number N.7128 and by CNPq of Brazil under contract number 216026/2013-0 T. Gokhool, R. Martins and P. Rives are with INRIA Sophia-Antipolis, France {tawsif.gokhool, renato-jose.martins, patrick.rives}@inria.fr Noëla Despré is with Airbus Defence and Space, Toulouse, France noela.despre@astrium.eads.net estimates are available. Eventually, to maximise the benefits of these two complementary approaches, it is needed to strike the right balance between metric and topological maps which leads us to the term topo-metric maps as in [START_REF] Badino | Visual Topometric Localization[END_REF].

In this context, accurate and compact 3D environment modelling and reconstruction has drawn increased interests within the vision and robotics community over the years as it is perceived as a vital tool for Visual SLAM techniques in realising tasks such as localisation, navigation, exploration and path planning [START_REF] Bylow | Real-Time camera tracking and 3D Reconstruction using signed distance functions[END_REF]. Planar representation of the environment achieve fast localisation [START_REF] Fernández-Moral | Fast place recognition with plane-based maps[END_REF], but they ignore other useful information in the scene. Space discretisation using volumetric methods such as signed distance functions [START_REF] Newcombe | KinectFusion: Real-Time Dense Surface Mapping and Tracking[END_REF], with combined sparse representations using octrees [START_REF] Steinbrücker | Large-scale Multi-Resolution Surface Reconstruction from RGB-D Sequences[END_REF], though they have received widespread attention recently due to their reconstruction quality do however present certain caveats. Storage capacity and computational burden restrict their application to limited scales. To provide alternatives, a different approach consisting of multi-view rendering was proposed in [START_REF] Merrell | Real-Time Visibility-Based Fusion of Depth Maps[END_REF]. Depth maps captured in a window of n views are rendered in a central view taken as the reference frame and are fused based on a stability analysis of the projected measurements. On a similar note, recent works undertaken in [START_REF] Meilland | On unifying keyframe and voxel based dense visual SLAM at large scales[END_REF] adopted a multi-keyframe fusion approach by forward warping and blending to create a novel synthesized depth map.

When exploring vast scale environments, many frames sharing redundant information clutter the memory space considerably. Furthermore, performing frame to frame registration introduces drift in the trajectory due to uncertainty in the estimated pose as pointed out in [START_REF] Kerl | Dense Visual SLAM for RGB-D Cameras[END_REF]. The idea to retain keyframes based on predefined criteria proves very useful within a sparse skeletal pose graph, which is the chosen model for our compact environment representation.

II. METHODOLOGY AND CONTRIBUTIONS

Our aim is concentrated around building ego-centric topometric maps represented as a graph of keyframes, spread by spherical RGBD nodes. A locally defined geo-referenced frame is taken as an initial model which is then refined over the course of the trajectory by neighbouring frames' rendering. This not only reduces data redundancy but also help in suppressing sensor noise whilst contributing significantly in drift reduction. A generic uncertainty propagation model is devised and leaned upon a data association framework for discrepancy detection between the measured and observed data. We build upon the above two concepts to introduce the notion of landmark stability over time. This is an important aspect for intelligent 3D points selection which serve as better potential candidates for subsequent inter frame to keyframe motion estimation task.

This work is directly related to two previous syntheses of [START_REF] Dryanovski | Fast visual odometry and mapping from rgb-d data[END_REF] and [START_REF] Meilland | On unifying keyframe and voxel based dense visual SLAM at large scales[END_REF]. The former differs from our approach in the sense that it is a sparse feature based technique and only consider depth information propagation. On the other hand, the latter aboard a similar dense based method to ours but without considering error models and a simpler Mahalanobis test defines the hypothesis for outlier rejection. Additionally, we build on a previous work of [START_REF] Meilland | A spherical robot-centered representation for urban navigation[END_REF] which constitutes of building a saliency map for each reference model, by adding the concept of stability of the underlying 3D structure.

Key contributions of this work are outlined as follows:

• development of a generic spherical uncertainty error propagation model, which can be easily adapted to other sensor models (e.g. perspective RGBD cameras) • a coherent dense outlier rejection and data fusion framework relying on the proposed uncertainty model, yielding more precise spherical keyframes • dynamic 3D points are tracked along the trajectory and are pruned out by skimming down a saliency map The rest of this paper is organised as follows: basic concepts are first introduced and shall serve as the backbone to allow a smooth transition between respective sections. Then the uncertainty error model is presented followed by the fusion stage. Subsequently, the problem of dynamic 3D points is tackled leading to a direct application of the saliency map. An experimental and results section verifies all of the above derived concepts before wrapping up with a final conclusion and some perspectives for future work.

III. PRELIMINARIES

An augmented spherical image S = {I, D} is composed of I ∈ [0, 1] m×n as pixel intensities and D ∈ R m×n as the depth information for each pixel in I. The basic environment representation consists of a set of spheres acquired over time together with a set of rigid transforms T ∈ SE(3) connecting adjacent spheres (e.g. T ij lies S j and S i ) -this representation is well described in [START_REF] Meilland | Dense visual mapping of large scale environments for real-time localisation[END_REF].

The spherical images are encoded in a 2D image and the mapping between the image pixel coordinates p and depth to cartesian coordinates is given by g : (u, v, 1) → q, g(p) = ρq S (p), with q S being the point representation in the unit spherical space S 2 and ρ = D(p) is radial depth. The inverse transform g -1 corresponds to the spherical projection model.

Point coordinates correspondences between spheres are given by the warping function w, under observability conditions at different viewpoints. Given a pixel coordinate p * , its coordinate p in another sphere related by a rigid transform T is given by a screw transform from the 3D q = g(p * ) followed by a spherical projection:

p = w(p * , T) = g -1 " [I 0]T -1 » g(p * ) 1 -« , ( 1
)
where I is a (3 × 3) identity matrix and 0 is a (3 × 1) zeros vector. In the following, spherical RGBD registration and keyframe based environment representations are introduced.

A. Spherical Registration and Keyframe Selection

The relative location between raw spheres is obtained using a visual spherical registration procedure [12] [7]. Given a first pose estimate T of T = TT(x), the sphere to sphere registration problem consists of incrementally estimating a 6 degrees of freedom (DOF) pose x ∈ R 6 , through the minimization of a robust hybrid photometric and geometric error cost function

F S = 1 2 e I 2 I + λ 2 2 e ρ 2 D
, which can be written explicitly as:

FS = 1 2 P p * W I (p * ) ' ' 'I " w(p * , b TT(x)) " -I * `p * ´' ' ' 2 + λ 2 2 P p * W D (p * ) ' ' 'n T " g(w(p * , b TT(x))) -b TT(x)g * (p * ) "' ' ' 2 (2)
where w(•) is the warping function as in (1), λ is a tuning parameter to effectively balance the two cost functions, W I and W D are weights relating the confidence of each measure and n is the normal computed from the cross product of adjacent points from q(p * ). The linearization of the over-determined system of equations ( 2) leads to a classic iterative Least Mean Squares (ILMS) solution. Furthermore for computational efficiency, one can choose a subset of more informative pixels (salient points) that yield enough constraints over the 6DOF, without compromising the accuracy of the pose estimate. This simple registration procedure applied for each sequential pair of spheres allows to represent the scene structure, but subjected to cumulative VO errors and scalability issues due to long-term frame accumulation. The idea to retain just some spherical keyframes, based on predefined criteria, proves very useful to overcome these issues and generate a more compact graph scene representation. A criteria based on differential entropy approach [START_REF] Kerl | Dense Visual SLAM for RGB-D Cameras[END_REF] has been applied in this work for keyframe selection.

IV. SPHERICAL UNCERTAINTY PROPAGATION AND MODEL FUSION

Our approach to topometric map building is an egocentric representation operating locally on sensor data. The concept of proximity used to combine information is evaluated Fig. 2. Full system pipeline mainly with the entropy similarity criteria after the registration procedure. Instead of performing a complete bundle adjustment along all parameters including poses and structure for the full set of close raw spheres S i to the related keyframe model S * , the procedure is done incrementally in two stages.

The concept is as follows: primarily, given a reference sphere S * and a candidate sphere S, the cost function in ( 2) is employed to extract T = TT(x) and the entropy criteria is applied for a similarity measure between the tuple {S * , S}. While this metric is below a predefined threshold, the keyframe model is refined in a second stage -warping S and carrying out a geometric and photometric fusion procedure is composed of three steps:

• warping S and its resulting model error propagation • data fusion with occlusions and outlier rejection • wise 3D point selection based on stable salient points which are detailed in the following subsections. The full system pipeline is depicted in figure [START_REF] Blanco | A tutorial on se(3) transformation parameterizations and on-manifold optimization[END_REF].

A. Warped Sphere Uncertainty

Warping the augmented sphere S generates a synthetic view of the scene S w = {I w , D w }, as it would appear from a new viewpoint. This section aims to represent the confidence of the elements in S w , which clearly depends on the combination of an apriori pixel position, the depth and the pose errors over a set of geometric and projective operations -the warping function as in [START_REF] Badino | Visual Topometric Localization[END_REF]. Starting with D w , the projected depth image is

Dw(p * ) = Dt(w(p * , T)) and Dt(p) = p qw(p, T) ⊤ qw(p, T) with qw(p, T) = " [I 0]T » g(p) 1 -« (3)
The uncertainty of the final warped depth σ 2 Dw then depends on two terms Σ w and σ 2

Dt ; the former relates to the error due to the warping w of pixel correspondences between two spheres and the latter, to the depth image representation in the reference coordinate system σ 2

Dt . Before introducing these two terms, let's represent the uncertainty due to the combination of pose T and a cartesian 3D point q errors. Taking a first order approximation of q = g(p) = ρq S , the error can be decomposed as:

Σq(p) = σ 2 ρ qSq ⊤ S + ρ 2 Σq S = σ 2 ρ ρ 2 g(p)g(p) ⊤ + ρ 2 Σ g(p)/ρ (4)
The basic error model for the raw depth is proportional to its fourth degree itself: σ 2 ρ ∝ ρ 4 , which can be applied to both stereopsis and active depth measurement systems (for instance see [START_REF] Khoshelham | Accuracy and resolution of kinect depth data for indoor mapping applications[END_REF] for details). The next step consists of combining the uncertain rigid transform T with the errors in q. Given the mean of the 6DOF x = {t x , t y , t z , θ, φ, ψ} in 3D+YPR form and its covariance Σ x , for q w (p, T) = Rq + t = Rg(p) + t , Σq w (p, T) = Jq(q, x)ΣqJq(q, x) ⊤ + JT(q, x)ΣxJT(q, x)

⊤ = RΣqR ⊤ + MΣxM ⊤ , (5)
with Σ q as in (4) and the general formula of M is given in [START_REF] Blanco | A tutorial on se(3) transformation parameterizations and on-manifold optimization[END_REF]. The first term Σ w using ( 5) and ( 1) is given by:

Σw(p * , T) = J g -1 (qw(•))Σq w (•)J g -1 (qw(•)) ⊤ (6)
where J g -1 is the jacobian of the spherical projection (the inverse of g) and (•) = (p * , T -1 ). The second term expression for the depth represented in the coordinate system of the reference sphere using the warped 3D point in ( 5) and ( 3) is straightforward

σ 2 Dt (p * , T) = JD t (qw(p * , T))Σq w (p * , T)JD t (qw(p * , T)) ⊤ (7)
with

J Dt (z) = (z ⊤ / √ z ⊤ z).

The uncertainty index σ 2

Dw is then the normalized covariance given by:

σ 2 Dw (p) = σ 2 Dt (p)/(q w (p, T) ⊤ q w (p, T)) 2 (8)
Finally, under the assumption of Lambertian surfaces, the photometric component is simply I w (p * ) = I(w(p * , T)) and it's uncertainty σ 2 I is set by a robust weighting function on the error using a Huber's M-estimator as in [START_REF] Meilland | Dense visual mapping of large scale environments for real-time localisation[END_REF].

B. Spherical Keyframe Fusion

Before combining the keyframe reference model S * with that of the transformed observation S w , a probabilistic test is performed to exclude outlier pixel measurements from S w , allowing fusion to occur only if the raw observation agrees with its corresponding value in S * .

Hence, the tuple A = {D * , D w } and B = {I * , I w } are defined as the sets of model predicted and measured depth and intensity values respectively. Finally, let a class c : D * (p) = D w (p) relate to the case where the measurement value agrees with its corresponding observation value. Inspired by the work of [START_REF] Murarka | Building Local Safety Maps for a Wheelchair Robot using Vision Lasers[END_REF], the Bayesian framework for data association leads us to:

p(c|A, B) = p(A,B|c)p(c) p(A,B) (9)
Applying independence rule between depth and visual properties and assuming a uniform prior on the class c (can also be learned from supervised techniques), the above expression simplifies to:

p(c|A, B) ∝ p(A|c)p(B|c) (10)
Treating each term independently, the first term of equation ( 10 Since equation ( 10) can be viewed as a likelihood function, it is easier to analytically work with its logarithm in order to extract a decision boundary. Plugging p(A|c) and p(B|c) into [START_REF] Khoshelham | Accuracy and resolution of kinect depth data for indoor mapping applications[END_REF] and taking its negative log gives the following decision rule for an inlier observation value:

(Dw(p) -D * (p)) 2 σ 2 Dw (p) + σ 2 D * (p) + (Iw(p) -I * (p)) 2 σ 2 Iw (p) + σ 2 I * (p) < λ 2 M , (12)
relating to the square of the Mahalanobis distance. The threshold λ 2 M is obtained by looking up the χ 2 2 table. Ultimately, we close up with a classic fusion stage, whereby depth and appearance based consistencies are coalesced to obtain an improved estimate of the keyframe sphere. Warped values that pass the test in [START_REF] Meilland | On unifying keyframe and voxel based dense visual SLAM at large scales[END_REF] are fused up by combining their respective uncertainties as follows:

I * k+1 (p) = W I k (p)I * k (p) + ΠI (p)Iw(p) W I k (p) + ΠI (p) , D * k+1 (p) = W D k (p)D * k (p) + ΠD(p)Dw(p) W D k (p) + ΠD(p) (13)
for the intensity and depth values respectively and weight update:

W I k+1 = W I k + ΠI and W D k+1 = W D k + ΠD (14)
where Π I (p) = 1/σ 2 Iw (p) and Π D (p) = 1/σ 2 Dw (p) from the uncertainty warp propagation in sec. IV-A.

C. Dynamic 3D points filtering

So far, the problem of data fusion of consistent estimates in a local model has been addressed. But to improve the performance of any model, another important aspect of any mapping system is to limit if not completely eliminate the negative effects of dynamic points. These points exhibit erratic behaviours along the trajectory and as a matter of fact, they are highly unstable. There are however different levels of "dynamicity" as mentioned in [START_REF] Konolige | Towards lifelong Visual Maps[END_REF]. Points/landmarks observed can exhibit a gradual degradation over time, while others may undergo a sudden brutal change -the case of an occlusion for example. The latter being considerably apparent in indoor environments where small viewpoint changes can trigger a large part of a scene to be occluded. Other cases are observations undergoing cyclic dynamics (doors opening and closing). Whilst the above-mentioned behaviours are learned in clusters [START_REF] Konolige | Towards lifelong Visual Maps[END_REF], in this work, points with periodic dynamics are simply evaluated as occlusion phenomena.

The probabilistic framework for data association developed in the section IV-B is a perfect fit to filter out inconsistent data. 3D points giving a positive response to test equation ( 12) are given a vote 1, or otherwise attributed a 0. This gives rise to a confidence map C * i (k) which is updated as follows:

C * i (k + 1) = ( C * i (k) + 1, if λ (95%) M < 5.991 0, otherwise (15)
Algorithm 1 3D Points Pruning using Saliency map

Require: ˘S * sal , C * i (k), N, n return Optimal Set A 10% ∈ S * sal Initialise new A for i=S * sal (p * ) = 1 to S * sal (p * ) = max do compute p (occur) (p * i ) compute γ k+1 (p * i ) compute p (stable) (p * i ) if p (stable) (p * i ) ≥ 0.8 then A[ i ] ←-p * i end if if length(A[ i ]) ≥ A 10% then break end if end for
Hence, the probability of occurence is given by:

p(occur) = C * i (k + N ) N , ( 16
)
where n is the total number of accumulated views between two consecutive keyframes. p(occur), though it gives an indication on how many times a point has been tracked along the trajectory, it can however not distinguish between noisy data or an occlusion. Treading on a similar technique to that adopted in [START_REF] Johns | Generative Methods for Long-Term Place Recognition in Dynamic Scenes[END_REF], a Markov observation independence is imposed. In the event that a landmark/3D point has been detected at time instant k, it is most probable to appear again at k + 1 irrespective of its past history. On the contrary, if it has not been re-observed, this may mean that the landmark is quite noisy/unstable or has been removed indeterminately from the environment and has little chance to appear again. These hypotheses are formally translated as follows:

γ k+1 (p * ) = ( 1, if p * k = 1 (1 -p(occur)) n , otherwise (17)
Finally, the overall stability of the point is given as:

p(stable) = γ k+1 p(occur) (18)

D. Application to Saliency map

Instead of naively dropping out points below a certain threshold, for e.g, p(stable) < 0.8, they are better pruned out of a saliency map [START_REF] Meilland | A spherical robot-centered representation for urban navigation[END_REF]. A saliency map, S * sal , is the outcome of careful selection of the most informative points, best representing a 6 degree of freedom pose, x ∈ R 6 , based on a Normal Flow Constraint spherical jacobian. The underlying algorithm is outlined in algorithm [START_REF] Badino | Visual Topometric Localization[END_REF]. The green and red sub-blocks in figure (3) represent the set of inliers b TT(x)

`( b TT(x) ´-1 set B set A A ∩ B w A -A 10% ∩ B w A 10%
Fig. 3. Saliency map skimming and outliers respectively, while the yellow one corresponds to the set of pixels which belong to the universal set {U : U = A ∪ B w } but which have not been pruned out. This happens when the Keyframe criteria based on an entropy ratio α [7] [START_REF] Kerl | Dense Visual SLAM for RGB-D Cameras[END_REF] is reached. The latter is an abstraction of uncertainty related to the pose x along the trajectory, whose behaviour shall be discussed in the results section.

The novelty of this approach compared to the initial work of [START_REF] Meilland | A spherical robot-centered representation for urban navigation[END_REF] is two-fold. Firstly, the notion of uncertainty is incorporated in spherical pixel tracking. Secondly, as new incoming frames are acquired, rasterised and fused, the information content of the initial model is enriched and hence the saliency map needs updating. This gives a newly ordered set of pixels to which is attributed a corresponding stability factor. Based on this information, an enhanced pixel selection is performed consisting of pixels with a greater chance of occurence in the subsequent frame. This set of pixel shall then be used for the forthcoming frame to keyframe motion estimation task. Eventually, between an updated model at time t 0 and the following re-initialised one, at t n , an optimal mix of information sharing happens between the two.

V. EXPERIMENTS AND RESULTS

A new sensor for a large field of view RGBD image acquisition has been used in this work. This device integrates 8 Asus Xtion Pro Live (Asus XPL) sensors as shown in figure (4)(left) and allows to build a spherical representation. The chosen configuration offers the advantage of creating full 360 • RGBD images of the scene isometrically, i.e. the same solid angle is assigned to each pixel. This permits to apply directly some operations, like point cloud reconstruction, photo consistency alignment or image subsampling. For more information about sensor calibration and spherical RGBD construction, the interested reader is referred to [START_REF] Fernández-Moral | Extrinsic calibration of a set of range cameras in 5 seconds without pattern[END_REF].

Our experimental test bench consists of the spherical sensor embarked on a mobile experimental platform and driven around in an indoor office building environment for a first learning phase whilst spherical RGBD data is acquired online and registered in a database. In this work, we do not address the problem of the real time aspect of autonomous navigation and mapping but rather investigate ways of building robust and compact environment representations by capturing the observability and dynamics of the latter. Along this streamline, reliable estimates coming from sensor data based on 3D geometrical structure are combined together to serve as a useful purpose for later navigation and localisation tasks which can thereafter be treated with better precision online.

Our algorithm has been extensively tested on a dataset 1 of around 2500 frames covering a total distance of around 60m. Figures 5(a), (b) illustrates the trajectories obtained from two experimented methods, namely; RGBD registration without (method 1) and with (method 2) keyframe fusion in order to identify the added value of the fusion stage. The visible discrepancy between the trajectories for the same pathway 1 video url:(http://www-sop.inria.fr/members/Tawsif.Gokhool/icra15.html) highlights the effects of erroneous pose estimation that occurs along the trajectories. This is even more emphasized by inspecting the 3D structure of the reconstructed environment as shown by figures (4)(centre), (4)(right) where the two images correspond to the sequence with and without fusion; method 1 and method 2 respectively. In detail, reconstruction with method 1 demonstrates the effects of duplicated structures (especially surrounding wall structures) which is explained by the fact that point clouds are not perfectly stitched together on revisited areas due to inherent presence of rotational drift, that is more pronounced than translational drift. However, these effects are well reduced by the fusion stage though not completely eliminated as illustrated in figure 4(right). The red dots on the reconstruction images are attributed to the reference spheres initialised along the trajectories using the keyframe criteria described in section III-A. 270 key spheres were recorded for method 1 while only 67 were registered for method 2.

Finally, comparing figures 5(a), (b), the gain in compactness for method 2 is clearly demonstrated by the sparse positioning of keyframes in the skeletal pose graph structure. Figure 5(c) depicts the behaviour of our keyframe selection entropy-based criteria α . The threshold for α is generally heuristically tuned. For method 1, a preset value of 0.96 was used based on the number of iterations to convergence of the cost function outlined in section III-A. With the fusion stage, the value of α was allowed to drop to 0.78 with generally faster convergence achieved. Figure 5(d) confirms the motion estimation quality of method 2 as it exhibits a lower error norm across frames as compared to method 1. Figure [START_REF] Badino | Visual Topometric Localization[END_REF] depicts the quality of the depth map before and after applying the filtering technique. The colours in the figure represent normal orientations with respect to the underlying surfaces.

VI. CONCLUSION

In this work, a framework for hybrid metric and topological maps in a single compact skeletal pose graph representation has been proposed. Two methods have been experimented and the importance of data fusion has been highlighted with the benefits of reducing data noise, redundancy, tracking drift as well as maintaining a compact environment representation using keyframes.

Our activities are centered around building robust and stable environment representations for lifelong navigation and map building. Future investigations will be directed towards gathering more semantic information about the environment to be able to further exploit the quality and stability of the encompassing 3D structure.

Fig. 1 .

 1 Fig. 1. Normal consistency between raw and improved spherical keyframe model of one of the nodes of the skeletal pose graph. The colours in the figure encode surface normal orientations, where better consistence is achieved when local planar regions depicts the same colour

  ) is devised as p(A|c) = p(D w (p)|D * (p), c), whereby marginalizing over the true depth value ρ i leads to: p(Dw(p)|D * (p), c) = Z p(Dw(p)|ρi, D * (p), c)p(ρi|D * (p), c)dρi (11) Approximating both probability density functions as Gaussians, the above integral reduces to: p(A|c) = N (D w (p)|D * (p), σ Dw (p), σ D * (p)). Following a similar treatment, p(B|c) = N (I w (p)|I * (p), σ Iw (p), σ I * (p)).

Fig. 4 .

 4 Fig. 4. Left:Multi RGBD acquisition system, centre and right:Reconstruction quality comparison (bird-eye view)

Fig. 5 .

 5 Fig.5. Performance comparison between method 1: RGBD spherical registration without keyframe fusion and method 2 : RGBD spherical registration with fusion. (plots in red and blue refer to methods 1 and 2 respectively)

ACKNOWLEDGMENT

The authors extend their warmest gratitude to Maxime Meilland and Panagiotis Papadakis for their critical insights and suggestions in improving the material of this paper.


Hrishikesh Deshpande

Pierre Maurel

Christian Barillot

ADAPTIVE DICTIONARY LEARNING FOR COMPETITIVE CLASSIFICATION OF MULTIPLE SCLEROSIS LESIONS

Keywords: Sparse Representations, Adaptive Dictionary Learning, Magnetic Resonance Imaging

Sparse representations allow modeling data using a few basis elements of an over-complete dictionary and have been used in many image processing applications. We propose to use a sparse representation and an adaptive dictionary learning paradigm to automatically classify Multiple Sclerosis (MS) lesions from MRI. In particular, we investigate the effects of learning dictionaries specific to the lesions and individual healthy brain tissues, which include White Matter (WM), Gray Matter (GM) and Cerebrospinal Fluid (CSF). The dictionary size plays a major role in data representation but it is an even more crucial element in the case of competitive classification. We present an approach that adapts the size of the dictionary for each class, depending on the complexity of the underlying data. The proposed algorithm is evaluated on clinical data demonstrating improved classification.

INTRODUCTION

An important variety of natural images can be represented as linear combinations of a few atoms of a learned dictionary. Over the last few years, researchers have shown that the sparse representation of the images in such a way produces promising results in image denoising and classification, including face recognition and texture classification [START_REF] Elad | Image denoising via learned dictionaries and sparse representation[END_REF][START_REF] Wright | Robust face recognition via sparse representation[END_REF][START_REF] Mairal | Supervised dictionary learning[END_REF]. Recently, its applications in disease detection have started evolving [START_REF] Weiss | Multiple sclerosis lesion segmentation using dictionary learning and sparse coding[END_REF][START_REF] Deshpande | Detection of multiple sclerosis lesions using sparse representations and dictionary learning[END_REF]. In this paper, we present a novel method using sparse representations and dictionary learning framework to classify MS lesions.

Multiple sclerosis is an autoimmune disease affecting the central nervous system and is characterized by the structural damages of axons. MRI is the best paraclinical method for the diagnosis of MS and treatment efficacy. Manual segmentation of MS lesions, however, is a laborious and time consuming task, and is prone to high intra-and inter-expert variability. Several MS lesion segmentation methods have been proposed over the last decades, with an objective of handling large variety of MR data and which can provide results that correlate well with expert analysis [START_REF] Lladó | Segmentation of multiple sclerosis lesions in brain MRI: A review of automated approaches[END_REF]. These methods either use lesion properties or tissue segmentation to help lesion segmentation.

We investigate the performance of automatic lesion classification by taking into account the tissue specific information. In the past, Weiss et al. [START_REF] Weiss | Multiple sclerosis lesion segmentation using dictionary learning and sparse coding[END_REF] treated lesions as outliers and achieved MS lesion segmentation using a single dictionary learned with the help of all image patches. Deshpande et al. [START_REF] Deshpande | Detection of multiple sclerosis lesions using sparse representations and dictionary learning[END_REF] reported improvements by learning class specific dictionaries for the lesion and healthy tissue patches. However, we can further enrich this model by using dictionaries specific to the different tissue types, such as WM, GM and CSF, as opposed to learning a single dictionary for healthy tissue patches. We explore the fact that various tissues as well as lesions appear in different intensity patterns in distinct MR modality images. For example, WM appears as the brightest tissue in T1-weighted image, but the darkest in T2-weighted images. Therefore, learning class specific dictionaries for individual tissues should further discriminate between lesion and non-lesion classes. We consider multiple MR modality images and learn the tissue specific dictionaries for WM, GM, CSF and lesion classes, to better classify MS lesions.

The dictionaries learned for each class are aimed at better representation of an individual class. However, if there exists differences in the data-complexity between classes, the relative under-or over-representation of either class will lead to worse classification. One idea for better classification could be to learn the dictionaries with adaptive sizes, in order to take into account the data variability between different classes. Thus, in addition to the dictionary learning strategy mentioned above, we also investigate the effect of modifying the dictionary sizes, leading to the proposition of adaptive dictionary learning. The basic idea is to learn the class specific dictionaries which are better adapted to the data and also complexity of the data. Previous works [START_REF] Deshpande | Detection of multiple sclerosis lesions using sparse representations and dictionary learning[END_REF][START_REF] Ramirez | An MDL framework for sparse coding and dictionary learning[END_REF] have also reported the effects of dictionary size in image classification. We investigate this in the particular case of classification.

We first describe sparse coding and dictionary learning in Section 2. The methodology is explained in Section 3, followed by results in Section 4 and conclusion in Section 5.

SPARSE CODING AND DICTIONARY LEARNING

Sparse coding is the process of finding a sparse coefficient vector a ∈ R K for representing a given signal x ∈ R N using a few atoms of an over-complete dictionary D ∈ R N ×K . It is CONFIDENTIAL. Limited circulation. For review only.

Preprint submitted to 2015 IEEE International Symposium on Biomedical Imaging.

Received November 24, 2014.

given by min a a 0 s.t. x -Da 2 2 ≤ ε, where . 0 denotes the l 0 norm and ε is the error in representation. Replacement of l 0 norm with the l 1 norm also results in sparse solution [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. The sparse coding problem can then be given by min

a x -Da 2 2 + λ a 1 , (1)
where λ balances the trade-off between the error and sparsity. For a set of signals {x i } i=1,.,m , we can find a dictionary D from the underlying data, such that each signal is represented by a sparse linear combination of its atoms. It is given by min

D,{ai} i=1,..,m m i=1 x i -Da i 2 2 + λ a i 1 .
(

) 2
The optimization is an iterative two-step process: Sparse coding with fixed D and the dictionary update with fixed a.

PROPOSED APPROACH

MR images from each patient are first processed for the removal of noise and non-brain tissues. The images are then registered and image patches of predefined size are extracted. These patches, after intensity normalization, are assigned to WM, GM, CSF or lesion class. For this purpose, we use manual lesion segmentation images and tissue segmentation obtained using SPM [START_REF] Ashburner | Unified segmentation[END_REF]. We then learn the dictionaries using training data for each class and perform patch based classification and voxel-wise classification for the given test subject.

Each step is described in the subsections below.

Pre-processing

The noise introduced during MR acquisition is removed using non-local means and intensity inhomogenity correction [START_REF] Tustison | N4ITK: Improved N3 bias correction[END_REF][START_REF] Coupe | Optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images[END_REF]. To ensure the spatial correspondence, the images are registered with respect to T1-w MPRAGE volume [START_REF] Wells | Multi-modal volume registration by maximization of mutual information[END_REF] and are processed further to extract the intra-cranial mask [START_REF] Smith | Fast robust automated brain extraction[END_REF]. This limits the further analysis to the brain region.

Patch Extraction and Labeling

We divide the whole intracranial MR volume for each patient into 3-D patches, with a patch around every 2 voxels in each direction. The individual image patches of each MR modality are then flattened, concatenated together and are normalized.

The patches are then labelled as either WM, GM, CSF or lesion class. If the number of lesion voxels in the corresponding image block of the manual segmentation image exceeds the pre-defined threshold T L , we assign this patch to the lesion class. Otherwise, the patch is assigned to either WM, GM or CSF class, depending on the maximum number of voxels that belong to the corresponding class.

The labelled image patches are then divided into training and test data, and the experiments are performed by following Leave-One-Subject-Out-Cross-Validation (LOSOCV).

Classification using Adaptive Dictionary Learning

The class specific dictionaries are learned using training data of the corresponding class [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF]. The subsections below describe different strategies adopted while learning these dictionaries and the scheme of patch based classification. In every method, we obtain the sparse codes for the test patches using Eq (1), knowing the dictionary D k for the class k [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF]. Similar to [START_REF] Deshpande | Detection of multiple sclerosis lesions using sparse representations and dictionary learning[END_REF], we use a single dictionary to represent the healthy class: WM, GM and CSF together. We learn the class specific dictionaries D k of same size, for the healthy (k = 1) and lesion (k = 2) class. For a given test patch y i , the classification is performed by calculating the sparse coefficients a k i for each class and the test patch is then assigned to class c such that c = argmin

k y i -D k a k i 2 2 .
(3)

Two-Class: Different Dictionary Size (2-C DDS)

In this method, we allow different dictionary sizes for the healthy and lesion classes, in order to take into consideration the data variability between classes and the number of training samples in each class. For accounting to more variability and larger training size, we allow larger dictionary size for the healthy class and obtain MS lesion classification.

Four-Class: Same Dictionary Size (4-C SDS)

The fact that every tissue, WM, GM and CSF, appears in different intensity pattern in each MR modality, using a single dictionary for representing tissues might not be as effective as learning separate dictionaries for each tissue. Each dictionary is representative of its own class and reconstruction of the test data using true class dictionary would give the minimum reconstruction error. We thus perform classification based on reconstruction error in a similar manner, as mentioned above.

Four-Class: Different Dictionary Size (4-C DDS)

Here, we experiment with different dictionary sizes for WM, GM, CSF and lesion classes, for similar reasons mentioned in Section 3.3.2.

Voxel-wise Classification

As already stated, we classify the patches centered around every 2 voxels in each direction. For voxel-wise classification, we assign each voxel to either of the classes by using majority voting. The voxel is assigned to a class using majority votes of all patches that contain the voxel.

In the context of lesion classification, we finally record the number of voxels that belong to True Positives (TP), False CONFIDENTIAL. Limited circulation. For review only.

Preprint submitted to 2015 IEEE International Symposium on Biomedical Imaging.

Received November 24, 2014.

Negatives (FN) or False Positives (FP), and calculate sensitivity (SEN)= T P T P +F N and Positive Predictive Value (PPV) = T P T P +F P .

EXPERIMENTS AND RESULTS

We evaluated the proposed approach on MR dataset of 14 MS patients acquired on a Verio 3T Siemens scanner. T1-w MPRAGE, T2-w, PD and FLAIR modalities were chosen for the analysis. The volume size for T1-w MPRAGE and FLAIR was 256×256×160 and voxel size was 1×1×1 mm 3 . For T2w and PD, the volume size was 256×256×44 and voxel size was 1×1×3 mm 3 . Annotations of the lesions were carried out by an expert radiologist. One subject with strong MR artifacts was excluded from the analysis.

For labeling patches, we used threshold T L = 6, as mentioned in Section 3.2. For each subject, the number of lesion patches varied from 1K to 30K, depending on the lesion load, whereas the average number of patches for WM, GM and CSF were 50K, 90K and 30K respectively. Using LOSOCV, the patch size of 5×5×5 and the sparsity parameter λ = 0.95 were found to be optimal choices. We used these parameters for all experiments described next.

The voxel-wise classification results for the methods described above are shown in Table 1. Using same dictionary size of 5000 for healthy and lesion class, as denoted by method (a), results in high sensitivity but very low PPV. Considering more variability associated with the healthy class, we then use different dictionary sizes, 5000 for healthy and 1000 for lesion class. This drastically reduces FP, improving PPV, but also decreases sensitivity. We further enrich this model by learning tissue specific dictionaries for WM, GM, CSF. Using 4-class approach and the same dictionary size of 2000 for each tissue and lesion class gives better compromise between sensitivity and PPV. Finally, 4-class method with different dictionary sizes, 2000 each for WM, GM and CSF, and 1000 for lesion class, increases both the mean sensitivity and PPV values, as compared to the 2-class method with different dictionary sizes. The mean dice-scores for methods (b) and (d), as referred to in Table 1, are 43.46% and 48.20%, respectively. Their comparison also shows a significant difference in PPV and dice-scores, with respective p-values of 0.0086 and 0.0298. This confirms that learning class specific dictionaries for each tissue improves the classification.

The classification performance for various dictionary sizes in the 4-class methods is shown in Table 2. The results for method (c), in row X, show that the sensitivity and PPV increase until the dictionary sizes are increased to 2000. The dictionaries capture more details with the increase in their size, but sensitivity reduces thereafter, possibly because of the over-fitting in either class dictionary. In rows Y and Z, we fix the size of either lesion or tissue dictionaries and vary the size of the other. It is observed that the best results, in terms of both sensitivity and PPV, are obtained for the dictionary Pat.

( size of 2000 for each tissue and 1000 for lesion class. It can also be observed that it is the relative dictionary size that drives the classification and is more important than just the absolute dictionary size for each class.

It is crucial to adapt the size of the dictionaries to better control the classification. For such purpose, we analyzed the data using Principal Component Analysis (PCA), which gives an estimate of the intrinsic dimensionality of the data. It was found that for each tissue, GM, WM and CSF, approximately twice as many eigenvalues are required for an arbitrary proportion of the cumulative data variance (90%, 95% or 98%), as that required for the lesion data. Thus, as exhibited by method (d), it supports our adaption of the dictionary size for each tissue twice that for the lesion dictionary. For method (b), although the experimentally observed optimal dictionary size ratio of 5 was not found with PCA, the factor of 2 indicated by PCA still favors using higher dictionary size for the healthy class. One reasoning behind this might be the inability of PCA to analyze the non-linearity in the data. The intrinsic dimensionality estimation for this highly non-linear data could be further point of investigation.

In Figure 1, we show the results for patient 6, for all the methods discussed above. It can be seen that the methods (a) and (c), each using the dictionaries with the same size, suffer with many FP. The over-detections are reduced in methods (b) and (d), respectively, which use the adapted dictionary sizes. However, the 2-class method (b) has many FN. Including tissue specific information in the approach and using dictionaries of the adapted sizes results in significant improvement in the lesion classification with reduction in both FP and FN. This is shown in method (d), supporting our claim that the Table 2. Effect of dictionary sizes in classification. First two columns indicate the dictionary size for each tissue and lesion, whereas the last two columns indicate the sensitivity and PPV. Row X: Method (c) with various dictionary sizes, Rows Y and Z: Method (d) with fixed dictionary size for either lesion or tissues, while varying dictionary sizes for the other class. method with the tissue specific dictionaries and adapted dictionary sizes is a better choice over the 2-class methods and those using the same dictionary size. 1. Classification image is superimposed on FLAIR MRI. TP are in red, FP are in cyan, FN are in green.

CONCLUSION

An approach for MS lesion classification based on dictionary learning and sparse representations is introduced this work. Learning tissue specific dictionaries, instead of using a single dictionary for the non-lesion class, has shown clear improvement in the lesion classification. We have also demonstrated the effectiveness of adapting the dictionary sizes for better amplification of the differences among classes, hence improv-ing the classification. If performing PCA on input data can successfully adapt the dictionary size for the classification, it is not as much efficient when the classes represent more a mixture of different tissues. Knowing the limitation of PCA to handle only linear data, a future work will be to use the intrinsic dimension estimation techniques, which can better analyze complexity of the non-linear data.

3. 3 . 1 .

 31 Two-Class: Same Dictionary Size (2-C SDS)

Fig. 1 .

 1 Fig. 1. Comparison of MS lesion classification methods, as referred in Table 1. Classification image is superimposed on FLAIR MRI. TP are in red, FP are in cyan, FN are in green.

Table 1 .

 1 Voxel-wise classification results using: (a) Two dictionaries with 5000 atoms for healthy and lesion classes (b) 5000 atoms for healthy and 1000 atoms for lesion class dictionary. (c) Four dictionaries with 2000 atoms for WM, GM, CSF and lesion class, (d) 2000 atoms for WM, GM and CSF, and 1000 atoms for lesion class dictionary.

	)

Preprint submitted to 2015 IEEE International Symposium on Biomedical Imaging.Received November 2014.

Preprint submitted to 2015 IEEE International Symposium on Biomedical Imaging.Received November 24, 2014.


Manikandan Bakthavatchalam

François Chaumette

Omar Tahri

An Improved Modelling Scheme for Photometric Moments with Inclusion of Spatial Weights for Visual Servoing with Partial Appearance/Disappearance

come   L'archive ouverte pluridisciplinaire

An Improved Modelling Scheme for Photometric Moments with Inclusion of Spatial Weights for Visual Servoing with Partial Appearance/Disappearance

Manikandan Bakthavatchalam, Franc ¸ois Chaumette, Omar Tahri Abstract-Visual servoing based on photometric data is of great interest since it does not necessitate any image processing or visual tracking steps. A vital issue in such methods is the change in the image resulting from the appearance and disappearance of portions of the scene from the camera field-ofview during the visual servo. In this paper, we propose a spatial weighting scheme to counter this problem. A general model for photometric moments that allows for inclusion of spatial weights is presented. A custom weighting function is then chosen and specific properties of its derivative are exploited in order to develop the interaction matrix in analytical form. The resultant effects on the invariance properties are discussed. Experiments have been performed to validate the visual servoing with the proposed weighting scheme on a 4 dof robot arm.

I. INTRODUCTION

Visual Servoing(VS) encompasses a whole spectrum of methods wherein the objective is to control the motions of a robot with a stream of images acquired from a camera sensor [START_REF] Hill | Real time control of a robot with a mobile camera[END_REF]. In Image-based Visual Servoing (IBVS), the error to be regulated (e = ss * ) is defined in the image space [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF] [START_REF] Iwatsuki | A new formulation of visual servoing based on cylindrical coordinate system[END_REF]. The s are the visual features built from image measurements from the projections of geometric (points, straight lines) or volumetric (spheres, cylinders) entities in the scene. In [START_REF] Silveira | Direct visual servoing: Vision-based estimation and control using only nonmetric information[END_REF] and [START_REF] Benhimane | Homography-based 2d visual tracking and servoing[END_REF], visual features derived from projective homography are used. This is estimated using a photo-geometric image registration algorithm, which by itself is a tracking problem. All these methods therefore require reliable visual tracking and/or matching processes.

Despite decades of research in VS, its potential for realworld applications is still not maximal. Among several factors that impede this expansion, dependence of VS on visual tracking procedures has been identified as a major factor [START_REF] Gridseth | Bringing visual servoing into real world applications[END_REF] [START_REF] Collewet | Visual servoing set free from image processing[END_REF]. Methods like [START_REF] Nayar | Subspace methods for robot vision[END_REF] used Principal Component Analysis (PCA) to transform image intensities to a lower dimensional representation. However, it was computationally prohibitive to expand them to several degrees of freedom [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF]. Pure photometric visual servoing [START_REF] Collewet | Photometric visual servoing[END_REF] directly utilised the luminance information in the image. It offers an excellent accuracy due to its usage of highly redundant features. It is also highly non-linear which reduces its convergence domain and produces non-intuitive robot trajectories.

Image moments are attractive for IBVS and shown to be better than geometric features because of their decoupling properties [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF] [START_REF] Tahri | Point-based and region-based image moments for visual servoing of planar objects[END_REF]. But their applications remain limited because of the restriction that well-segmented regions of interest or a discrete point set should always be available in the image. [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF] tried to remove this restriction with photometric image moments defined from the image plane intensities. While no visual tracking or matching were required, all the decoupling properties of moments were retained. However, the modelling did not take into account the resultant effects induced by portions of the image which are not common between the current image and the next (including the image learnt from the desired robot pose). Since extraneous regions can disturb the stability and convergence of the VS, it is important to solve this vital issue. Interestingly, this issue arising from the uncommon image portions was pointed out earlier in [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF].

Image moments [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF] [13] use polynomial kernels and are a specific case of KBVS [START_REF] Kallem | Kernelbased visual servoing[END_REF]. Kernels of finite support were used in [START_REF] Swensen | Empirical characterization of convergence properties for kernel-based visual servoing[END_REF] perhaps to deal with this issue. The nonavailability of analytical models for the interaction matrix in KBVS (and in learning-based approaches like [START_REF] Nayar | Subspace methods for robot vision[END_REF]) make conducting further analyses impossible.

Works like [START_REF] García-Aracil | Continuous visual servoing despite the changes of visibility in image features[END_REF] and [START_REF] Durand | Management of visual signal loss during image based visual servoing[END_REF] deal with feature visibility and signal loss in IBVS, problems of somewhat similar nature. But image points are used as features and so not applicable to dense photometric methods. To the best of our knowledge, solutions to this problem do not currently exist.

The main idea in this paper is to introduce spatial weights in image moments so that the areas where the appearance and disappearance occur (at the periphery) have as less influence as possible when compared to the central parts of the image during the visual servoing.

Our paper is organized as follows: Sec II presents the modelling aspects to obtain the interaction matrix of weighted photometric moments and the selection of the spatial weighting function. Sec III discusses the properties for the most commonly used visual features. Finally, validation results are presented in Sec IV.

II. MODELLING OF WEIGHTED PHOTOMETRIC MOMENTS

We define a weighted photometric moment (WPM) of order (p + q) as follows

m pq = π x p y q w (x) I (x, t) dx dy (1)
where x = (x, y) is a spatial point where the intensity I(x, t) is measured at time t, π denotes the image support and w(x) is a spatial weighting function.

A. General case

To maintain generality, we first develop an analytical form of the interaction matrix for an abstract spatial weighting function w (x) and then for a specific family of functions in Sec II-B. In the following derivation, we follow the same line of reasoning as in [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF], only emphasizing the changes introduced by the weighting function, wherever appropriate. We recall that ṁpq = L mpq v c where L mpq is the interaction matrix of the moments defined in (1) and v c = (v x , v y , v z , ω x , ω y , ω z ) is the camera velocity screw. By taking the derivative of (1), we have

ṁpq = π x p y q w (x) İ (x, t) dx dy (2)
where İ (x, t) = L I v c , L I being the interaction matrix of the luminance feature [START_REF] Collewet | Photometric visual servoing[END_REF]. Let us note that this result is obtained from the optic flow constraint equation (OFCE) with the brightness constancy hypothesis [START_REF] Horn | Determining optical flow[END_REF]. After following the same developments as in [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF], we obtain

                             L vx mpq = A m ∇x p+1,q + B m ∇x p,q+1 + C m ∇x p,q L vy mpq = A m ∇y p+1,q + B m ∇y p,q+1 + C m ∇y p,q L vz mpq = -A m ∇x p+2,q -B m ∇x p+1,q+1 -C m ∇x p+1,q -A m ∇y p+1,q+1 -B m ∇y p,q+2 -C m ∇y p,q+1 L ωx mpq = -m ∇x p+1,q+1 -m ∇y p,q -m ∇y p,q+2
L ωy mpq = m ∇x p,q + m ∇x p+2,q + m ∇y p+1,q+1

L ωz mpq = -m ∇x p,q+1 + m ∇y p+1,q (3)
If we let Q = x p y q w(x) I(x) and P = 0 in (5), we have

π x p y q w(x) ∂I ∂x + x p y q ∂w ∂x I(x) (6)
+p x p-1 y q w(x) I(x) dx dy = ∂π x p y q w(x) I(x)dy So, equation (4a) becomes

m ∇x pq = - π
x p y q ∂w ∂x I(x) + p x p-1 y q w(x)I(x) dx dy + ∂π

x p y q w(x, y)I(x, y)dy

(7)
Using the definition of the weighted moments in (1), we subsequently have

m ∇x pq = -p m p-1,q - π x p y q ∂w ∂x I(x) dx dy + ∂π x p y q w(x, y)I(x, y)dy (8)
The last term has to be evaluated over the image border ∂π.

In [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF], this term was removed with the assumption that the image borders are of zero intensity. This hypothesis is invalid when there are dynamic image regions which enter and leave during the visual servo. Here, we define the weighting function so that w(x, y) = 0 ∀(x, y) ∈ ∂π and so, the last term vanishes again but without the zero border assumption. Further, if we take P = x p y q w(x, y)I(x, y) and Q = 0 in Green's theorem and make the same series of developments as before, m ∇y pq can be obtained. The terms in (4) related to the interaction matrix of the WPM become

m ∇x pq = -p m p-1,q - π x p y q ∂w ∂x I(x, y)dxdy (9a)
m ∇y pq = -q m p,q-1π

x p y q ∂w ∂y I(x, y)dxdy (9b)

Substituting [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF] in Eq(3), we see that the interaction matrix is dependent not only on weighted moments m pq , but also has additional terms composed of π x p y q ∂w ∂x I (x). For this, obviously the first order spatial derivatives should exist for w(x) 1 . Next, we select an appropriate weighting function such that a closed-form expression for the interaction matrix is obtained as functions of the WPM.

B. Weighting scheme

The weighting function should attribute maximal importance to the pixels in the area around the image centre, smoothly reducing it outwards towards the image edges. As pointed out earlier, the first order spatial derivatives of w(x) should exist ∀(x, y) ∈ π. The standard logistic function l(x) = 1 1+exp -x was initially considered for the weighting scheme. With this choice however, it was not possible to obtain the interaction matrix in an analytical form as functions of the weighted photometric moments. That is why in this paper, we propose to use functions with the general structure given in [START_REF] Collewet | Photometric visual servoing[END_REF], where the exponential is raised to the negative power of a polynomial function.

F(x) = exp -p(x) with (10) p(x) = a 0 + a 1 x + 1 2 a 2 x 2 + 1 3 a 3 x 3 + ... + 1 n a n x n
Functions of this structure possess the interesting property that their derivatives can be expressed in terms of the function itself. The derivative of F(x) is given by

F (x) = exp -p(x) p (x) with (11) p (x) = a 1 + a 2 x + a 3 x 2 + ... + a n x n-1
This in turn allows to obtain the interaction matrix as a function of the weighted moments (see Eq15). We propose the following weighting function with the form as in ( 10)

w(x, y) = K exp -a(x 2 +y 2 ) 2 ( 12
)
where K is the maximum value of the weighting function and a can be used to vary the areas which receive maximum and 0 weights respectively (See Fig. 1). The polynomial function (x2 + y 2 ) 2 to which the exponential is raised preserves the invariance property of the original image moments (i.e., moments without weights) to rotations around the optic axis. This can be confirmed from L ωz mpq in [START_REF] Swensen | Empirical characterization of convergence properties for kernel-based visual servoing[END_REF]. We now continue

∂w ∂x = -4ax(x 2 + y 2 ) w(x) (13a) ∂w ∂y = -4ay(x 2 + y 2 ) w(x) (13b)
Substituting ( 13) into ( 9), we obtain

m ∇x pq = -p m p-1,q + 4a (m p+3,q + m p+1,q+2 ) (14a)
m ∇y pq = -q m p,q-1 + 4a (m p,q+3 + m p+2,q+1 ) (14b)
By insertion of ( 14) into (3), the interaction matrix of photometric moments L mpq weighted with the radial function in [START_REF] Tahri | Point-based and region-based image moments for visual servoing of planar objects[END_REF] can be obtained:

L mpq = L vx mpq L vy mpq L vz mpq L ωx mpq L ωy mpq L ωz mpq (15)
with:

L vx mpq = A (-(p + 1) m p,q + 4 a (m p+4,q + m p+2,q+2 )) + B (-pm p-1,q+1 + 4 a (m p+3,q+1 + m p+1,q+3 ))
+ C (-pm p-1,q + 4 a (m p+3,q + m p+1,q+2 ))

L vy mpq = A (-qm p+1,q-1 + 4 a (m p+3,q+1 + m p+1,q+3 )) + B (-(q + 1) m p,q + 4 a (m p,q+4 + m p+2,q+2 )) + C (-qm p,q-1 + 4 a (m p,q+3 + m p+2,q+1
))

L vz mpq = A[(p + q + 3) m p+1,q -4 a (m p+5,q + 2m p+3,q+2 + m p+1,q+4 )] + B[(p + q + 3) m p,q+1 -4 a (m p+4,q+1 + 2m p+2,q+3 + m p,q+5 )] + C[(p + q + 2) m p,q -4 a (m p+4,q + 2m p+2,q+2 + m p,q+4 )] L ωx mpq = (p + q + 3)m p,q+1 + qm p,q-1 -4 a(m p+4,q+1 + 2 m p+2,q+3 + m p,q+3 + m p+2,q+1 + m p,q+5 ) L ωy mpq = -pm p-1,q -(p + q + 3)m p+1,q + 4 a(m p+3,q + m p+1,q+2 + m p+5,q + 2 m p+3,q+2 + m p+1,q+4 ) L ωz mpq = pm p-1,q+1 -qm p+1,q-1
This interaction matrix is expressed as a matrix sum

L mpq = L 1 + 4aL 2 ( 16
)
where L 1 has the same form as the non-weighted moments in [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF] (but the moments are now weighted as in (1)). L 2 is tied directly to the weighting function. To compute L mpq , moments m tu of order t + u <= p + q + 5 are required. So, the choice of the weighting function dictates the maximum moment order required for computing the interaction matrix.

In the case of a uniform weighting function, the above interaction matrix reduces 2 exactly to the case in [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF]. In L mpq , it is to be noted that L ωz mpq has no new term contributed by the weighting function. Thus, the invariance of the weighted moments to optic axis rotations is retained. This outcome was of course desired from the structure of the weighting function in [START_REF] Tahri | Point-based and region-based image moments for visual servoing of planar objects[END_REF]. Considering the other components in L mpq , we observe additional terms induced by the weighting function. As a result, invariants to scale and 2D translations built from the classical moments are not invariant anymore when using the WPM. This loss of invariance is an issue for selecting visual features to control the full 6 dof of a system.

III. ANALYSIS OF THE LOWER ORDER WEIGHTED MOMENTS

A. Weighted photometric area moment

With WPM, the zeroth order moment can be interpreted as the weighted photometric area. Its interaction matrix can be derived from [START_REF] Swensen | Empirical characterization of convergence properties for kernel-based visual servoing[END_REF] simply by taking p = q = 0 and by using K = 4a for brevity. We obtain

                 L vx m00 = -m 0,0 A + K [Aα 1 + Bα 2 + Cλ 1 ] L vy m00 = -m 0,0 B + K [Aα 2 + Bα 3 + Cλ 2 ] L vz m00 = m 00 (3/Z g -C) -K [Aλ 3 + Bλ 4 + C(α 1 + α 3 )] L ωx m00 = 3m 00 y g -K (λ 2 + λ 5 ) L ωy m00 = -3m 00 x g + K (λ 1 + λ 3 ) L ωz m00 = 0 (17)
with

           α 1 = m 4,0 + m 2,2 λ 1 = m 3,0 + m 1,2 α 2 = m 3,1 + m 1,3 α 3 = m 0,4 + m 2,2 λ 2 = m 0,3 + m 2,1 λ 3 = m 1,4 + 2m 3,2 + m 5,0 λ 4 = m 0,5 + 2m 2,3 + m 4,1
For the original non-weighted moments, the area m 00 is a very interesting visual feature mainly related to the translation motions in z. It is invariant to (i) rotations around the optic axis and (ii) to planar translations when the camera and target planes are parallel (A = B = 0). With WPM, from the terms in [START_REF] Durand | Management of visual signal loss during image based visual servoing[END_REF], we indeed observe that the invariance to optic axis rotations is preserved thanks to our choice of the weighting function. The influence on planar translations can be studied by taking A = B = 0 in the corresponding terms of the interaction matrix. Then, we will have L vx m00 = K C(m 3,0 + m 1,2 ) and L vy m00 = K C(m 0,3 + m 2,1 ). Clearly, third order moments appear, which means that the invariance property is generally lost. From the experimental results however, this loss did not have any strong effect on the convergence of the servo. On the contrary, it is to be noted that if these terms are neglected, it will introduce a gross error into the modelling and the visual servoing might not converge at all.

B. Weighted Centre of gravity

The ratio between the first and zeroth order moments gives the centre of gravity coordinates in the image x g = m 10 /m 00 and y g = m 01 /m 00 . The interaction matrix of x g can be obtained analytically using

L xg = L m1,0 -x g L m0,0 /m 0,0 .                                      L vx xg = -1/Z g + K m 0,0 A 1 + B 2 + C(α 1 -x g λ 1 ) L vy xg = K m 0,0 A 2 + B 1 + C(α 2 -x g λ 2 ) L vz xg = x g /Z g + 4 (An 20 + Bn 11 ) - K m 0,0 A(φ 1 -x g λ 3 ) + B(φ 2 -x g λ 4 ) + C(λ 5 -x g (α 1 + α 3 ))                          L ωx xg = (x g y g + 4n 1,1 ) + K m 0,0 (φ 2 + α 2 ) -x g (λ 2 + λ 4 ) L ωy xg = -(1 + x 2 g + 4n 20 ) + K m 0,0 (φ 1 + α 1 ) -x g (λ 1 + λ 5 ) L ωz xg = y g (18)
with

               i = δ i -x g α i δ 1 = m 5,0 + m 3,2 δ 2 = m 4,1 + m 2,3 φ 1 = m 2,4 + 2m 4,2 + m 6,0 φ 2 = m 5,1 + 2m 3,3 + m 1,5 λ 5 = m 5,0 + 2m 3,2 + m 1,4
and n i,j = µ i,j /m 0,0 where µ i,j = m i,jm 0,0 x i g , y j g . From the second term L vy xg in the interaction matrix, we can notice that translational velocities in y will modify the value of x g even when the camera and target planes are parallel. Hence x g is not invariant to motions in y. This loss of invariance to translations was unfortunately expected from the spatial weighting. We can notice the weighting function did not alter the behaviour of the centre of gravity w.r.t optic axis rotations. In a similar manner, the interaction matrix of the other coordinate y g can also be obtained. But since it does not present any particularly interesting characteristics, it is omitted here.

C. Photometric orientation feature

The orientation α is an interesting feature that can be calculated from the centred second order moments. α = 1 2 arctan 2µ 11 µ 20µ 02 [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF] where µ pq are the weighted centred moments. The exact analytical form of L α is not given here due to space constraints. It is straight-forward to compute however as shown in [11] [12]. As with the case of the centre of gravity and area, the behaviour with respect to optic axis rotations is not modified. That is, the last term of the interaction matrix is equal to -1, indicating that it is linearly related to rotational motion around the optic axis. This linear relation shows α is well-suited to control motion along this dof.

IV. VALIDATION RESULTS

We present in this section both simulation and experimental results that validate the proposed approach. Software for this research was built on top of the ViSP [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF] library. A calibrated camera is mounted on the end-effector in an eye-in-hand configuration and a motionless, planar target is considered. For simulations as well as the experiments with the robot, a SCARA-type (3T+1R) actuation with only the 3D translations and rotation around the optical axis was considered. This is the same subset of dof considered in KBVS [START_REF] Kallem | Kernelbased visual servoing[END_REF]. The partial velocity screw v c = (v x , v y , v z , ω z ) is computed using the classical IBVS control law [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF]:

v c = -λ L -1 s (s -s * ) ( 20
)
with s = (x n , y n , a n , α), where

(x n = x g • a n , y n = y g • a n )
is the normalized weighted centre of gravity, a n = Z * m * 00 m00

is the normalized area moment, (x g , y g ) the centre of gravity and α is the orientation as defined in [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]. These features are used in several works including [START_REF] Tahri | Point-based and region-based image moments for visual servoing of planar objects[END_REF] and [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF] (but from moments without weighting). x n , y n , a n are normalized versions of the features described in Sec III and λ is the gain for the control law. The weighting function ( 12) with K = 1 and a = 1000 was satisfactory (no precise tuning was required). Further, the spatial weights are computed only once in a single offline step and kept in a convenient matrix form.

A. Simulation Results

For the simulations, the desired pose is chosen 1.0m vertically above the target such that the camera and target planes are parallel. A gain of λ = 2.0 was used for faster convergence (not affected by lower values). The initial pose is chosen so as to obtain the image shown in Fig. 2b. From Figs. 2a and2b, we can observe that there is no distinguishable background and the ship structure remains obstructed in the initial image. The displacements required for convergence are a translation of c * t c = [2.41cm, 6.65cm, 10cm] and a rotation of R z (25 • ) around the optic axis.

1) Using L s = L s * : The control law in (20) was used with the interaction matrix computed once at the desired pose L s (s = s * ). Hence, the control action is linear in (ss * ) since L s * is constant and is simpler than in [START_REF] Collewet | Photometric visual servoing[END_REF]. The interaction matrix computed at the desired pose has a condition number of 2365. Discussion : From Eq.( 21), several inferences can be made. The feature a n is invariant to optic axis rotations and α has retained its property with respect to optic axis rotations (-1 in last column). But x n (resp. y n ) is no more invariant to translations in y (resp. x). This invariance loss was discussed in Sec.II-B. Also, a n and α in (21) are more reactive to errors in degrees of freedom different from the ones they are intended to control (translations along and around the optic axis resp). This consequence is due to the weighting function, which implies new features would have to be determined for a more optimal design. At this juncture, we would like to stress the importance of explicit analytical modelling which facilitates such analyses. Despite the loss in invariance, the visual servoing converged satisfactorily, especially with respect to the translation motions. The rotation error increases briefly before the final convergence. This is due to the extraneous regions appearing in the image during the initial stages of the servo. The 3D camera trajectory is not close to a straight line (See Fig. 2f) but an overall satisfactory behaviour was obtained. The camera traversed a distance of 14.11cm in this experiment. Comparison to non-weighted moments : With the same desired and initial poses, the visual servoing with nonweighted photometric moments proposed in [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF] was tested. As seen from Fig. 2d, the decrease in the errors for the nonweighted moments is not satisfactory at all and the servoing diverged. After about 40 iterations, the algorithm using the proposed WPM is driven near to convergence whereas for the non-weighted moments in [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF], there is a large error as seen from the difference image in Fig. 2g. The experiment was stopped when the error norm increased rapidly, otherwise leading to a jitter. It is evident that the control law in this case cannot handle the extraneous regions.

2) Using L s = L s : The same initial and desired poses as in Sec.IV-A.1 are considered. For the interaction matrix, it is computed at the current pose L s (s) for use in the control law (20). As seen from the results in Fig. 3a, a perfect The departure from a straight line behaviour is more in this case (see Fig. 3c) than the previous one (see Fig. 2f). This exponential error decrease and more deviation of the camera trajectory were both expected. This is because usage of the current interaction matrix ensures a nice exponential regulation of the errors but does not ensure the same for the Cartesian trajectory [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF].

B. Experimental Results

For the experiments, the Viper850 platform was configured for a SCARA-type actuation. The desired camera pose is chosen at approximately 40cm in front of the target. The displacement required between the initial and desired poses is c * t c = [-2.32cm, 0.97cm, -8.3cm] in translation and R z (-11 • ) rotation around the optic axis. Like in the simulations, these poses are chosen such that there is no distinction between what is the foreground and what constitutes the background in the initial and desired images. As for the control law, a gain of λ = 0.3 and the interaction matrix L s * at the desired configuration is used.

1) Experiment I: The image learnt from the desired pose is shown in Figure 4a. The image acquired from the initial pose is shown in Figure 4b. Clearly, the actress who is present in the left portion of the initial image is not present at all in the learnt image. This is a voluntary choice as explained previously. The interaction matrix at the desired position had a condition number of 818.39 with same invariance properties discussed in IV-A. The extraneous image regions have induced large variations in the error in a n during the initial iterations as seen from Fig. 4c. Nevertheless, the system is driven to convergence by the generated control velocities shown in Fig. 4d. The final error in the achieved robot pose is [-0.29mm, -2.21mm, -0.96mm] in translation and 0.1 • around the optic axis respectively. 2) Experiment II: In this experiment, the same initial and desired poses as Experiment IV-B.1 is chosen, but with a different photometric target. The images acquired at the desired and initial poses are shown in Figures 5a and5b respectively. The form of the interaction matrix is similar to that in (21). Its condition number is 4435.2, higher than in the previous experiment. Indeed, the exact values in the interaction matrix depends on the spatial distribution of the intensities in the image. Further, from Figure 5c, we find that the error in a n is more noisy than in the previous experiment. This is due to the fact that the noise effect depends on the area Experiment II : Visual Servoing on the Viper850 robot with weighted photometric moments under the influence of the weighting function. This can also be attributed in part to the uncontrolled lighting conditions which can contribute to noise in the pixels. The accuracy in this experiment is [0.84mm, -1.72mm, -0.95mm] in translations which is acceptable and an error of only 0.07 • around the optic axis. The proposed scheme was thus validated by several successful experiments despite the loss in invariance. We recall that our scheme does not require any image matching or visual tracking processes but only the computation of a set of weighted moments on the image.

The hypothesis about the information being present in the central image parts might not always hold. If the shared area between the current and learnt image is very low, the method might fail (and this is true for all VS methods in general). This approach is also affected by illumination variations, which is indeed a limitation.

V. CONCLUSIONS AND FUTURE WORK

In this paper, an enhanced modelling scheme for photometric moments that accommodates spatial weights was presented. This improved model took into account scenario of extraneous image regions. A carefully selected weighting function was introduced into the photometric moments formulation to obtain the interaction matrix in analytical form. Simulations and experiments with a real robot configured for SCARA-type actuation validate the proposed method. The weighted photometric moments preserved their behaviour with respect to rotations around the optic axis while behaviour with respect to other degrees of freedom were not retained. This affects the invariance properties of moment invariants, on which the features that control rotations around the image axes [START_REF] Tahri | Point-based and region-based image moments for visual servoing of planar objects[END_REF] are based. Establishment of this lost invariance is therefore an important problem and requires a considerable advance in theory, which is left to future work.

Fig. 1 .

 1 Fig. 1. Surface of the weighting function in (12) with K = 1, a = 650 with our developments with (12) as the weighting function. Its spatial derivatives are as follows:

Fig. 2 .

 2 Fig. 2. Servo Results for Sec.IV-A.1 (a) Desired image, (b) Initial image, Visual feature errors for (c) Weighted Photometric Moments and for (d) non-weighted moments [13], both shown in Multi-scale y-axis (xn, yn, an on left and α on right) (e) Camera velocities and (e) Camera 3-D trajectory for Weighted Photometric Moments, Difference image after 39 iterations (f) using non-weighted moments and (g) using Weighted Photometric Moments

Fig. 3 .

 3 Fig. 3. Visual Servoing Results with the proposed Weighted Photometric Moments using b Ls = Ls exponential convergence of the errors is achieved in this case. This concretely validates the correctness of all the modelling steps proposed in this paper. The camera traversed a distance of about 18.32cm, greater than the distance of 14.11cm traversed in the earlier experiment in Sec.IV-A.1.The departure from a straight line behaviour is more in this case (see Fig.3c) than the previous one (see Fig.2f). This exponential error decrease and more deviation of the camera trajectory were both expected. This is because usage of the current interaction matrix ensures a nice exponential regulation of the errors but does not ensure the same for the Cartesian trajectory[START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF].

Fig. 4 .

 4 Fig. 4. Experiment I : Visual Servoing on the Viper850 robot with weighted photometric moments

  Fig. 5.Experiment II : Visual Servoing on the Viper850 robot with weighted photometric moments

Note that if w(x, y) = 1 ∀(x, y), the gradients of w would vanish. Both m ∇x pq and m ∇y pq would then reduce to a single term (which was the case in[START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF] as well). That is, they reduce exactly to the case of photometric moments with no weight

This can be easily verified since w(x, y) = 1 would imply a = 0 and therefore Lm pq = L 1 and only moments of order (p + q + 1) are sufficient to compute it.