Guillaume Laville

Kamel Mazouzi

Christophe Lang

Nicolas Marilleau

Laurent Philippe

Using GPU for Multi-agent Multi-scale Simulations

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Scientists pay more and more attention to models. They aim at simulating real complex systems to understand them, for example to prevent phenomena such as disasters. Models thus (such as Agent Based Models -ABM-) become more and more descriptive: where scientists used to built up conceptual models (KISS model), they now create well described models (KIDS Model) [START_REF] Edmonds | from kiss to kids: An" anti-simplistic" modeling approach[END_REF]. Tomorrow the complex system community may tend to "world model", a fully descriptive and generic model that scientists can customize according to scientific questions [START_REF] Amouroux | KIMONO: using the modelling process as an aid for research orientation[END_REF]. Making these modeling evolutions possible requires a raise of computing power as the model executions must return a result in a reasonable delay.

If a simulation can be parallelized then increasing the computing power becomes a question of cost: we need to buy more cores to increase the size of the simulation. While standard processors (CPU) are still expensive, modern Graphical Processing Units (GPU) provide good execution performance for a lower cost. These GPUs also give the possibility to execute non-graphic programs using languages such as OpenCL or CUDA.

Typical simulators are however based on a sequential design and only use one core of the main processor. So taking advantage of the computer architecture implies to develop distributed simulators. In this context, Multi-Agent System (MAS) is an Nicolas Marilleau Institut de Recherche pour le Développement (IRD), France, e-mail: nicolas.marilleau@ ird.fr Christophe Lang • Guillaume Laville • Kamel Mazouzi • Laurent Philippe Institut FEMTO, CNRS / Université de Franche-Comté, France, e-mail: name.surname@ univ-fcomte.fr 1 interesting way to create this kind of models and simulators. ABM are indeed often used to simulate natural or collective phenomenons whose actors are too numerous or various to provide a unified algorithm describing the system evolution.

We propose in this paper to delegate part of the ABM execution to the graphical unit of the computer, based on our experience in parallelizing and implementing part of an ABM on a GPU. We then run it on both a standard CPU processor and on a GPU. We get very good performance results and advocate for a more wider use of the GPU in ABM in particular for simple agents. In Section 2 of the paper we present the work related to agent and parallelization. We give an overview of the sworm simulator in Section 3 and we detail its GPU implementation in 4. We present the experiments and results in Section 5 then we conclude on the possible generalization of our work.

Related works

To develop ABM, many frameworks are now available (e.g. Repast [START_REF] North | A declarative model assembly infrastructure for verification and validation[END_REF] or NetLogo [START_REF] Sklar | Netlogo, a multi-agent simulation environment[END_REF]). Only a few of them introduces distribution in agent simulation (Madkit [START_REF] Gutknecht | Madkit: a generic multi-agent platform[END_REF] or MASON [START_REF] Sean | Mason: A multi-agent simulation environment[END_REF]). In this context, parallel implementations are often based on threads using shared memory or on the integration of cluster libraries such as MPI. Even as multi-core or multi-socket setups become more and more readily available, these solutions stay limited to a relatively small number of parallel tasks, if not using a fully-fledged computing cluster.

Parallelizing a simulation is however complex as space and time constraints must be enforced. Time constraints are linked to the simulation execution. It is usually based on a synchronous execution of time steps by the agents and the environment. Distributing the simulation or delegating part of this execution to others processors as GPUs [START_REF] Bleiweiss | Multi agent navigation on the gpu[END_REF] must thus be carefully done to enforce the synchronism [START_REF] Marilleau | An agent based framework for urban mobility simulation[END_REF]. Space constraints are linked to the environment distribution. Using an environment on a set of computers leads to classical parallelism issues: data coherency if the environment is shared amongst all the computers or data exchanges if the environment is distributed. In the particular case of multi-scale simulations such as the Sworm simulation [START_REF] Blanchart | Sworm: an agent-based model to simulate the effect of earthworms on soil structure[END_REF] the environment may be used at different levels. This characteristic, especially in a fractal model, could be the key of the distribution. For instance, each branch of a fractal environment could be identified as an independent area and parallelized. In addition Fractal is a famous approach to describe multi-scale environment (such as soil) and its organization [START_REF] Bird | The psf and soil density scaling[END_REF].

Sworm

Sworm model aims at simulating soil functioning especially the role of macro-fauna on microbial activity and on carbon production. A first version presented in [START_REF] Blanchart | Sworm: an agent-based model to simulate the effect of earthworms on soil structure[END_REF] modified by these earthworms, are considered by biologist as hot spots of a microbial activity (an important process that transforms organic matters into carbon and so gives the fertility to the soil). The MIOR -MIcro-ORganisms-model tackles reproducing the microbial activity of these hot spots. MIOR focuses on a nearby 0.002 mm sized cube of soil while Sworm models 20 cm sized volume of soil.

Since these different models do not act at the same level of detail, this simulation is split into multiple, recursive levels corresponding to different scales. It uses a fractal pattern which allows lazy environment allocation: the different kinds of agents do not work at the same scale, depending of their respective sizes. As an example, worms tend to interact on a macroscopic scale (few millimeters), whereas microbial colonies work at a microscopic scale (10 -4 millimeters). Figure 1 illustrates the fractal data representation of Sworm. White cells (in the left part of the cell) represent soil cavities, black cells represent mineral matter and gray ones represent composite cells that must be further decomposed..

Scale n

Scale n + 1

Fig. 1 Sworm environment representation

A MIOR simulation is based on two types of agents. MMs -Meta-MIOR-represent colonies of carbon-consuming microbes naturally found in soil. OMs -Organic Matters-are deposits that represent the base resources randomly distributed in the modeled environment. Two main processes account for the model evolution, repre-senting the MM metabolism: breathing (conversion of mineral carbon to CO 2 ) and growth (fixation of carbon into microbial colonies). This second process only occurs if the colony needs are covered by the environment, as described in Algorithm 1.

Even with the sparse fractal representation, simulating the microscopic-scale evolution quickly becomes time-consuming on a traditional CPU unit, where much execution resources are already required for the macroscopic simulation (nearly one week per simulation on a personal computer).

GPU Implementation

CUDA and OpenCL are the main programming languages for GPU scientific programming. They are based on a C-like syntax with additions to support GPUspecific features. Various open-source bindings, such as JOCL [START_REF] Jocl | Java bindings for OpenCL[END_REF] and JCUDA [START_REF] Jcuda | Java bindings for CUDA[END_REF] exploit GPUs with the Java language. Since our implementation must be integrated in a Java-based implementation and must be portable JOCL is a natural choice.

Sworm requires two species of agents, cognitive ones (representing worms) and reactive ones (MIOR microbial colonies). The two kinds are not as good candidates for a GPU implementation. Reactive agents tend to rely on data parallelism with few, if any, kind of branching in the code to be executed. On the opposite, cognitive agents have a variety of possible behaviors at each iteration.

Our work therefore mainly focuses on the adaptation of the MIOR part of the Sworm simulation to a GPU architecture. Since the environment modeling and evolution if often a costly aspect in many ABM, theses changes can provide a tangible benefit on the simulator executions. By using simultaneously the GPU for the environment execution and the CPU for the cognitive agents we also improve the quality of modeling (replacement of the precomputed static carbon evolution graph by dynamically computed ABMs in our case).

OpenCL representation

Java can only be used to control the model initialization and the OpenCL executions. So a mapping of the Java Object-Oriented concepts to the C-like data structures of OpenCL is required to access and manipulate them in the GPU code. Here are some guidelines that can be used in this mapping process. The data aspect (type, attributes) can be represented as OpenCL structures. Methods can be implemented as device-local functions, and each phase of the simulation can have its own highlevel function calling them. Consecutive function calls with the same dimensions and subjects may be gathered into the same kernel (set of threads) to minimize the number of context switches between host and GPU.

By applying these mapping rules, four main data structures appears in the MIOR simulator: one array of MM structures for microbial colonies, one array of OM structures for organic deposits, one neighboring matrix used to store accessibility between MM and OM and one single environment used to store global parameters such as agent numbers and output data such as the quantity of CO 2 .

The design of the computing topology is of importance for the parallelism of the execution. Due to their simple behavior the reactive agents usually have limited interactions and thus allow a high level of parallelism. In our case, the action radius of each MM agent is determined by a fixed parameter of the environment, RA. All organic matter deposits within this distance should be accessible for the metabolism.

Since there are no dependencies between each distance calculation, we can use a two dimension kernel in OpenCL where each thread represents an (OM, MM) couple. This approach can be generalized with n dimension kernels in OpenCL depending on the modelized environment and on the GPU characteristics (usually one to three dimensions).

Data Representation and Data Dependencies

The limited interactions between reactive agents also impact the data representation as the density of interaction matrices may be low. In a typical MIOR simulation each agent is linked to a small (less than 10) number of OM on a total of some hundredth. The topology computation thus produces a low-density matrix. During the following evolution process, going through these matrices may induce numerous costly global memory accesses without actual use of the information. This issue can be addressed by creating dense continuous representations of the matrices. So on the topology stage we create two dense continuous representations of the neighboring matrix. The additional memory cost of these two data structures (one for the OM, one for the MM) can thereafter be balanced thanks to the model living steps.

Data dependency is a much difficult issue as it usually depends on the ABM characteristics. Some guidelines may be given but good results can only be obtained thanks to experience. Most of the time, as in our case, an algorithm adaptation will be needed. We present two cases of algorithm adaptation in this section.

In the MIOR simulation the data dependencies comes from the need to ensure fair access to carbon deposits for each microbial colony and thus to synchronize these accesses across the simulation. Since the colonies are randomly placed at the model initialization, no obvious geographic static locality can be extracted to split the carbon deposits in local memories for the sets of microbial colonies (MM). Therefore much data has to remain into the global slow GPU memory. A first implementation based on mutual execution and global memory synchronization thus resulted in GPU performances two order of scale slower than a sequential Java simulation. This clearly indicated the need for an algorithm adaptation. Our solution is to reduce the number of synchronization required for the execution of a given step of simulation. For that we split the living cycle process into three steps: 1. Scattering: carbon deposit resource is evenly scattered in parts across all microbial colonies accessing to it, 2. Living: each microbial colony consumes carbon for its breathing/growing process and produces CO 2 , 3. Gathering: each carbon deposit content is recomputed from the remaining carbon in each part.

There are also synchronization issues in the breathing and growth part of the initial algorithm. Since multiple MM can share a OM carbon deposit, two synchronization issues must be addressed to avoid introducing a bias in the simulation. First, a synchronization on the resources is needed as multiple MM agents should not modify the OM's properties at the same time. Second, the fair access to resources from the MM must be enforced and no OM should be over-exploited.

Data transfer between Host and GPU

Since the reference algorithm is sequentially implemented, a synchronization of the simulation steps must be implemented to prevent some agents to evolve faster than their peers. A first naive implementation may attempt to copy all modified data from and to Java objects between each simulation step. Data copy between the CPU system and the GPU-dedicated memory however uses a connexion limited to a fraction of the main memory bandwidth and with a much higher latency. So an obvious course of action to improve this is to allow the simulation to run several steps at once, or up to the termination avoiding these costly transfers. The only exception must be for values needed to detect simulation termination or convergence. The better improvements are achieved in the case where only the final state of the system is needed.

Proposed implementations

We propose three consecutive GPU implementations of the MIOR living process with an increased adaptation level. The first implementation, further referenced as GPU v1.0, is based on a straightforward adaptation of the CPU algorithm where the relations between agents are stored as a simple two-dimensions sparse matrix located in global memory. Each simulation step iterates over the whole matrix (including empty cells). The second GPU implementation, called GPU v2.0, replaces this matrix by a compressed one using the same representation format as described in [START_REF] Gómez-Luna | Parallelization of a video segmentation algorithm on cuda-enabled graphics processing units[END_REF]: neighbors are stored in a contiguous way in the matrix along with a precomputed count of neighbor agents. This representation reduces the required number of iterations and of memory accesses. To minimize global memory access overhead, the third GPU implementation, referenced as GPU v3.0, uses the GPU private memory to store often-used data such as carbon parts and MM neighboring information.

Experiments

To assess the performance of our work we compare the sequential version with our MIOR implementations on two platforms, representing what a researcher could except as dedicated computing hardware or personal computer graphic card. The first platform is a GPU node with two Intel Xeon X5550 CPU at 2.67GHz and a Tesla C1060 GPU card running at 1.3GHz (240 cores organized in 30 streaming multiprocessors). The second platform is a personal computer with an Intel Q9300 CPU at 2.5GHz and a mainstream GPU card: a GeForce 8800GT at 1.5GHz (112 cores organized in 14 streaming multiprocessors).

Figures 2 and3 give the average execution time of 50 simulations. The problem size is given by the scaling factor. A scaling factor of 6 means that the number of agents is multiplied by 6. At scale 1, the model contains 38 MM and 310 OM. We can note that small sized problems result in similar execution times across all the implementations. The GPU simulations do not have enough agents to benefit from the algorithm's parallelism. From scale 5 GPU versions v2.0 and v3.0 are clearly faster than the CPU or the simple GPU implementations. At scale 10, a ratio of 10 can be observed between the v3.0 GPU implementation and the CPU one.

Several important remarks must be reported. First, the optimized GPU versions give much better performance than the non-optimized GPU one. So it really matter to work on the algorithm parallelization. Second, even a mainstream desktop graphic card provides the same order of performance gain that a much more costly solution, the Tesla card. This is rather interesting as a much larger public could benefit from GPU adaptations. Third, the same code is run on the GeForce card as on the Tesla card without any modification. So the performance does not depend on the GPU characteristics and the same benefit could be expected on other graphic cards.

Conclusion

In this paper we have described an adaptation of an existing ABM simulation using GPU hardware. The first result of this work is that adapting the algorithm to a GPU architecture is possible for ABMs. This adaptation may provide a significant perfor-mance improvement without so much effort, or at least not much than for a standard multi-threaded parallelization. We have shown that this approach especially suits the case of multi-scale ABMs. We also have generalized this work to show the main issues to be addressed when parallelizing an ABM and how to take benefit from the GPU architecture. Last, we show that using a mainstream card, as the GPU card of a standard computer, can even lead to a significant performance improvement and avoid the use of a costly parallel cluster.

We are currently working on defining more general guidelines for adapting SMA simulations to the GPU architecture. The aim is to use both CPU and GPU at the same time to support multi-scale ABMs. With cognitive agents being run on the CPU and reactive agents or environment being run on the GPU, we could run larger simulations thanks to the performance improvement provided by this specific architecture. In this perspective, we are working on the design of a generic library to support reactive agent simulation on graphic cards.

Fig. 2 Fig. 3

 23 Fig.2CPU and GPU executions on a Tesla C1060

Acknowledgements Computations have been performed on the supercomputer facilities of the Mésocentre de calcul de Franche-Comté.


Aurelien Yol

Eric Marchand

Francois Chaumette

Keyvan Kanani

Thomas Chabot

VISION-BASED NAVIGATION IN LOW EARTH ORBIT

This paper presents and compares two vision-based navigation methods for tracking space debris in a low Earth orbit environment. This work is part of the RemoveDEBRIS 1 project. The proposed approaches rely on a frame to frame model-based tracking in order to obtain the complete 3D pose of the camera with respect to the target. The proposed algorithms robustly combine points of interest and edge features, as well as color-based features if needed. Experimental results are presented demonstrating the robustness of the approaches on synthetic image sequences simulating a CubeSat satellite orbiting the Earth. Finally, both methods are compared in order to specify what could be the best choice for the future of the RemoveDEBRIS mission.

INTRODUCTION

Since the beginning of the space era in 1957, the number of missions has incredibly increased. Unfortunately, it has also progressively generated a huge amount of debris that are still orbiting the Earth. Today, space agencies feel more and more concerned about Earth orbital environment cleaning since they realize those debris represent a big threat due to the high risk of collision that could, among safety issues, damage operational satellites [START_REF] Bonnal | Active debris removal: Recent progress and current trends[END_REF]. To investigate solutions to this issue, various projects aim to actively remove debris objects by catching them. Among these projects, is the RemoveDEBRIS project co-funded by the European Union and members of European space industry. In a first stage, the RemoveDEBRIS mission will consist of catching two miniaturized satellites (CubeSats DS-1 and DS-2 produced by the Surrey Space Centre) that will be previously ejected from a microsatellite called RemoveSat, itself released from the International Space Station (ISS). Then, the CubeSats will be used as targets instead of actual debris. This debris removal demonstration will

1
This work is supported by the European Commission FP7-SPACE-2013-1 (project 607099) « RemoveDEBRIS -A Low Cost Active Debris Removal Demonstration Mission », a consortium partnership project consisting of: Surrey Space Centre (University of Surrey), SSTL, Airbus DS, Airbus SAS, Airbus Ltd, ISIS, CSEM, INRIA, and Stellen-Bosh University.

include net capture, harpoon capture and vision-based navigation using standard camera and LiDAR [START_REF] Forshaw | An in-orbit active debris removal mission -REMOVEDEBRIS: Pre-Launch update[END_REF].

In the literature, several works propose approaches based on active sensors, such as LiDAR, for space rendezvous and docking missions [START_REF] Howard | Advanced video guidance sensor (AVGS) development testing, Defense and Security[END_REF]. Those approaches have proved to be efficient in situations where lighting conditions can be tricky due to directional sunlight resulting in high specularities and shadows. However, LiDAR are known to be generally expensive and rather heavy sensors. Furthermore, it also has a limited range and accuracy, and it cannot take the advantage of using potential texture information available on the target surface in order to improve the navigation process.

In this paper, we focus on the vision-based navigation part of the RemoveDEBRIS mission using a single standard RGB camera. Based on the knowledge of the 3D model of the target, estimating the complete 3D pose of the camera with respect to this target has always been an ongoing issue in computer vision and robotics applications [START_REF] Marchand | Pose estimation for augmented reality: a hands-on survey[END_REF]. For instance, regarding space applications, [START_REF] Oumer | Vision-based localization for onorbit servicing of a partially cooperative satellite[END_REF][START_REF] Petit | Tracking complex targets for space rendezvous and debris removal applications[END_REF] use model-based tracking approaches for space rendezvous with space target or debris. Common approaches try to solve this problem by using texture [START_REF] Bleser | Realtime 3d camera tracking for industrial augmented reality applications[END_REF], edge features [START_REF] Oumer | Vision-based localization for onorbit servicing of a partially cooperative satellite[END_REF][START_REF] Comport | Realtime markerless tracking for augmented reality: the virtual servoing framework[END_REF][START_REF] Drummond | Real-time visual tracking of complex structures[END_REF][START_REF] Petit | Tracking complex targets for space rendezvous and debris removal applications[END_REF], or color or intensity features [START_REF] Panin | An efficient and robust real-time contour tracking system[END_REF][START_REF] Prisacariu | Pwp3d: Real-time segmentation and tracking of 3d objects[END_REF][START_REF] Panin | Robust contourbased object tracking integrating color and edge likelihoods[END_REF]. The algorithms proposed in this paper provides a robust approach relying on a frame to frame tracking that align the projection of the 3D model of the target with observations made in the image by combining edges, point of interest and color-based features [START_REF] Pressigout | Real-time hybrid tracking using edge and texture information[END_REF][START_REF] Petit | Combining complementary edge, point and color cues in model-based tracking for highly dynamic scenes[END_REF]. Moreover, one of the presented method relies on the use of a 3D rendering engine to manage the projection of the model and to determine the visible and prominent edges from the rendered scene, while the other method is based on a hand-made 3D model but does not require any complex rendering process.

The paper is organized as follow: in a first stage, the general issues of the model-based tracking problem are recalled. Then, we describe how to combine edges features with keypoints and color-based features by considering two different model-based tracking approaches. The first approach only relies on the use of a CPU, while the second takes the advantage of using a graphics processing unit (GPU). Finally, experimental results are presented on synthetic image sequences where comparisons are done between both approaches.

MODEL-BASED TRACKING

Pose Estimation

As stated before, our problem is restricted to modelbased tracking where the considered 3D model is composed of a CAD representation of the target. In this paper, the model-based tracking part relies on a non-linear minimization of the reprojection error [START_REF] Marchand | Pose estimation for augmented reality: a hands-on survey[END_REF] where the general issue is to estimate the complete 3D pose of the camera with respect to the target 𝐌 # $ by minimizing the error ∆ between a set of measured data 𝑠 ' * (usually the position of a set of geometrical features in the current image) and the same features projected in the image plane according to the current pose 𝑠 ' 𝒓 :

∆(𝐫) = 𝜌(𝑠 ' 𝐫 -𝑠 ' * ) 0 1 234 ( 1
)
where 𝜌 is a robust estimator [START_REF] Huber | Robust Statistics[END_REF], which reduces the sensitivity to outliers (M-estimation) and 𝐫 is a vectorbased representation of

𝐌 # $ ( 𝐌 # $
is the homogeneous 4x4 matrix that represents the position of the object in the camera frame).

A common way of applying M-estimation is the Iteratively Reweighted Least Squares (IRLS) method. It converts the M-estimation problem into a leastsquares one where the error to be regulated to zero is defined as:

𝐞 = 𝐃(𝐬 𝐫 -𝐬 * ) ( 2
)
where 𝐃 is a diagonal weighting matrix whose each coefficient gives the confidence on each feature. Their computation is based on M-estimators [START_REF] Comport | Realtime markerless tracking for augmented reality: the virtual servoing framework[END_REF][START_REF] Drummond | Real-time visual tracking of complex structures[END_REF].

Then, by considering an iterative least square approach given by:

𝒗 = -𝜆 𝐃𝐋 𝒔 < 𝐞 ( 3
)
where 𝑣 is the velocity screw applied to the virtual camera position and 𝐋 𝒔 is the image Jacobian related to 𝐬 defined such that 𝐬 = 𝐋 𝒔 𝒗. The 3D pose of the camera can be iteratively updated from the computed displacement 𝒗 by applying an exponential map. More precisely, the resulting pose is given by:

𝐌 # >?@ $ = 𝐌 # > $ exp (𝒗) (4)
where 𝑘 denotes the number of iterations in the minimization process.

While considering model-based tracking, it is important to decide what are the features 𝐬 to be used in order to define the error criterion [START_REF] Oumer | Vision-based localization for onorbit servicing of a partially cooperative satellite[END_REF]. In this paper, we consider two model-based tracking methods that use a combination of geometrical edge-based features with keypoints features, and also color features along silhouette edges.

CPU-Based Model-Based Tracker

The first considered model-based tracking method relies only on a central processing unit (CPU), which is particularly appealing for space applications, thanks to the use of a simple hand-made CAD model. Furthermore, only edge-based features and keypoint features are used to define the error criterion [START_REF] Oumer | Vision-based localization for onorbit servicing of a partially cooperative satellite[END_REF].

Considering different kind of features allows to benefit from their complementarity and to overcome the limitations of a single feature-based approach. In order to use them in the non-linear process presented in the previous section, the error ∆ to be minimized has to be rewritten according to the considered features. For this first method, it is defined as:

∆ = 𝑤 F ∆ F + 𝑤 H ∆ H ( 5
)
where ∆ F refers to the geometrical edge-based error and ∆ H to the error relying on the keypoint features.

Respectively, 𝑤 F and 𝑤 H are the weighting parameters used to control the influence of each set of features.

Finally, to include the combination of the different features in the minimization process, the idea (as proposed in [START_REF] Petit | Combining complementary edge, point and color cues in model-based tracking for highly dynamic scenes[END_REF]) is to stack ∆ F and ∆ H in a single global error vector 𝐞 . Then, their corresponding Jacobian is also stacked and both are used as defined in ( 2) and (3).

Edge-based features

From the knowledge on the 3D model of the target, a potential option is to use its silhouette to define features to rely on. Thus and as defined in ( 5), the first features we consider for this method are edge-based features. As in [START_REF] Comport | Realtime markerless tracking for augmented reality: the virtual servoing framework[END_REF][START_REF] Drummond | Real-time visual tracking of complex structures[END_REF], the edges that are considered correspond to the projection of the CAD model of the target, using the current camera pose 𝐫 . More precisely, the edges corresponding to the visible polygons of the model are considered in the minimization process. That is why this approach requires a simplified version of the 3D model that only contains polygons with visible edges (in term of gradient in the image). Furthermore, the visibility process has to be well defined in order to not consider tracking is performed via a 1D search along the normal to the projection of the line 𝑙 ' 𝐫 which is associated to the point 𝐱 ' F (see Figure 1). As a result, the error criterion ∆ F referred as the edge-based features error is defined as:

∆ F (𝐫) = 𝜌 F (𝑑 P (𝑙 ' 𝐫 , 𝐱 ' F * )) 0 1 N 234 (6)
where 𝜌 F is a robust estimator and 𝑑 P (𝑙 ' 𝐫 , 𝐱 ' F * ) is the distance between the forward projection of the line 𝑙 ' 𝐫 of the 3D model of the target and the desired point features 𝐱 ' F * . Finally, the Jacobian of the distance of a point to a line is used (its analytic form is given in [START_REF] Comport | Realtime markerless tracking for augmented reality: the virtual servoing framework[END_REF]).

Keypoint-based features

a vision sensor and by considering that the target is potentially textured in the image, we have the possibility to use this information in order to improve the robustness of our navigation process. This is the reason why the considered methods also rely on keypoint-based features. Chosen keypoints are Harris corners that are tracked over the image sequence with the KLT algorithm [START_REF] Shi | Good features to track[END_REF]. Thus, the criterion error ∆ H refers here to the distance between 𝐱 ' H 𝐫 , the position of the initial extracted keypoints 𝐱 ' H with respect to the camera pose 𝐫, and their tracked position for each frame of the sequence 𝐱 ' H * (see Figure 1).

∆ H (𝐫) = 𝜌 H ( 𝐱 ' H 𝐫 -𝐱 ' H * ) 0 1 R 234 (7)
where 𝜌 H is once again a robust estimator.

Regarding 𝐱 ' H 𝐫 , let us recall that for this first method we consider a simplified version of the 3D model of the target where keypoint are supposed to be detected on planar polygons belonging to the model. Thus, the displacement between the initial feature position and the tracked one can be specified by an homography. Then 𝐱 ' H 𝐫 can be related to the initial position 𝐱 ' H of the keypoint by:

𝐱 ' H 𝐫 = α 𝐊 U V4 𝐇 𝐊 U 𝐱 ' H (8)
where 𝐊 U is the calibration matrix associated to the intrinsic camera parameters 𝜉 , and 𝐇 is the homography (up to a scale factor α) obtained from the 3D camera motion 𝐌 # YZY[ # > computed between the initial camera pose 𝐌 # YZY[ $ and the current camera pose 𝐌 # > $ :

𝐇 = 𝐑 ']'^ _ + 𝐭 ']'^ _ 𝑑 ']'^ 𝒏 ']'^b (9)
where Finally, the Jacobian considering the coordinates of points is used and given by:

𝐋 H = 𝐊 U - 1 𝑍 0 𝑥 𝑍 𝑥𝑦 -(1 + 𝑥 0 ) 𝑦 0 - 1 𝑍 𝑦 𝑍 (1 + 𝑦 0 ) 𝑥𝑦 -𝑥 (10)
where (𝑥, 𝑦) are the normalized coordinates of the image point 𝐱 ' H 𝐫 , and Z its corresponding depth obtained by:

1 𝑍 = 𝑑 -𝐭 ']'^ _ b 𝒏 ']'^ ']'( 𝐑 ']'^ _ 𝒏 ']'^) b 𝐱 ' H 𝐫 . (11)

GPU-Based Model-Based Tracker

For the second considered model-based tracking method, we take the advantage of using a low-end graphics processing unit (GPU) to improve the visionbased navigation. Furthermore, in addition to edgebased features and keypoint features that were used in the previous method, we also consider color-based features. As a result, the error ∆ to be minimized has to be rewritten as:

∆ = 𝑤 F ∆ F + 𝑤 H ∆ H + 𝑤 # ∆ # (12)
where ∆ F and ∆ H still refer to the geometrical edgebased and keypoint-based errors, with 𝑤 F and 𝑤 H their corresponding weighting parameters. ∆ # is the colorbased features error with its associated weighting parameter 𝑤 # .

Edge-based features

With respect to the previous method, the error criterion ∆ F remains the same (see [START_REF] Panin | An efficient and robust real-time contour tracking system[END_REF]). However, the determination of the edges to be considered is slightly different. In the approach proposed in [START_REF] Petit | Combining complementary edge, point and color cues in model-based tracking for highly dynamic scenes[END_REF] and used here, the CAD model of the target is rendered, according to the camera pose 𝐫 , by using a 3D rendering engine (in our case OpenSceneGraph). This allows to process the corresponding depth map of the scene in order to extract, via a Laplacian filter, the prominent edges of the projected model that will define the new set of initial point features { 𝐱 ' F } '34 M N .

Tracking phase remains the same but as we are now dealing with a much complex 3D model, points are not associated anymore with any polygon. As a result, control points are processed independently and a 3D line 𝑙 ' is computed for each 𝐱 ' F (in blue on Figure 2).

Keypoint-based features

As for the keypoints, the feature-based error criterion ∆ H remains the same (see [START_REF] Prisacariu | Pwp3d: Real-time segmentation and tracking of 3d objects[END_REF]). However, with respect to the previous method where polygons could contain several point features, we here have to take into account that every point may belong to a different polygon, due the potential complexity of the 3D model. Consequently, the way how 𝐱 ' H 𝐫 and its corresponding depth Z are computed has to be changed in order to avoid unnecessary time processing. As in [START_REF] Petit | Combining complementary edge, point and color cues in model-based tracking for highly dynamic scenes[END_REF] we take the advantage of using the GPU and more particularly the rendered depth map to back-project the initial detected keypoints { 𝐱 ' H } '34 M R on the model, which gives the set of corresponding 3D points { 𝐗 ' H } '34 M R . These 3D points are then projected with the value of the camera pose 𝐫 and each point 𝐱 ' H 𝐫 can be redefined as:

𝐱 ' H 𝐫 = 𝑝𝑟(𝐗 ' H , 𝐫) (13)
where 𝑝𝑟(𝐗 ' H , 𝐫) is the perspective projection operator.

Finally, the current depth value 𝑍 _ is computed by considering the current camera displacement

𝐌 # YZY[ # > . 𝐱 ' H _ = 𝐌 # YZY[ # > 𝐱 ' H ']'^ . ( 14
)
Note that the current depth value 𝑍 _ could also be obtained by rendering the scene at the current pose 𝐌 # > $ and then using the corresponding depth map. However, this solution has not been considered since it requires too much time to process the data.

Color-based features

One last improvement that has been added to this later approach in order to robustify the navigation process is to also consider color-based features. The idea is to use the prominent edges that have been extracted to compute the edge-based error criterion ∆ F in order to characterize the separation between the silhouette of the projected 3D model and the background by relying on color information [START_REF] Panin | Robust contourbased object tracking integrating color and edge likelihoods[END_REF][START_REF] Petit | Combining complementary edge, point and color cues in model-based tracking for highly dynamic scenes[END_REF].

The principle is to compute local color statistics (RGB means 𝐈 ' m and 𝐈 ' n and covariances 𝐑 ' m and 𝐑 ' n where O is the object and B the background) on both sides along the normal 𝒏 ' to the projected model silhouette edges { 𝐱 ' F } '34 M N , regularly sampled in 𝐱 ',q # points up to a distance 𝐿 (see Figure 2). These statistics are then mixed according to a fuzzy membership rule [START_REF] Petit | Robust visual detection and tracking of complex object: applications to space autonomous rendez-vous and proximity operations[END_REF], giving the mean 𝐈 ' (𝐫) and covariance 𝐑 ' (𝐫) where 𝐈 ' (𝐫) is defined as the measured color value for 𝐱 ',q # . Consequently, we can define the error 𝐞 ',q # 𝐫 as:

𝐞 ',q # 𝐫 = 𝐈 ' 𝐫 -𝐈(𝐱 ',q # ) (15)
and the general error criterion ∆ # can be re-written as:

∆ # (𝐫) = 𝜌 # ( 𝐞 ',q # 𝐫 b 𝐑 ' (𝐫) V4 𝐞 ',q # 𝐫 ) q ' [START_REF] Marchand | Pose estimation for augmented reality: a hands-on survey[END_REF] where 𝜌 # , as for 𝜌 H and 𝜌 F , is a robust estimator based on a Tukey M-Estimator.

For more accuracy and temporal smoothness, [START_REF] Panin | Robust contourbased object tracking integrating color and edge likelihoods[END_REF] and [START_REF] Petit | Combining complementary edge, point and color cues in model-based tracking for highly dynamic scenes[END_REF] also propose to introduce a temporal consistency by integrating the color statistics computed on the previous frame 𝐈 s for the silhouette edge point 𝐱 ' F 𝐫 _ at the first iteration of the minimization process. With a weighting factor ( 0 < α < 1 ), it gives 𝐞 ',q # 𝐫 = α𝐈 ' 𝐫 + (1 -𝛼)( 𝐈 ' 𝐫 s -𝐈 𝐱 ',q # ). For the detailed computation of the corresponding interaction matrix, see [START_REF] Petit | Robust visual detection and tracking of complex object: applications to space autonomous rendez-vous and proximity operations[END_REF].

Finally, let us note that in case an RGB camera would not be available, it is still possible to use this approach with intensity only images. Indeed, color statistics defined above could still be computed by considering one image component only, that is, the intensity level for grey images.

EXPERIMENTAL RESULTS

The

Structure of the target

The structure of the CubeSat DS-2 that will be used for the RemoveDEBRIS mission can be seen on Figure 3. It has avionics throughout the structure

{ I i O R i O I O , } L Measured position s i * R i B I i B ,
and four deployable panels at the bottom in the shape of a cross. The panels have no specific function except to make the CubeSat looking more like an actual satellite, so that the vision-based navigation is performed close to real conditions. Its dimensions (panels excluded) are 0.1m x 0.1m x 0.227m.

Case 1: case of a geosynchronous orbit

In the first sequence (see Figures 5 and6), the considered target is on a geosynchronous orbit. It has an angular speed that creates rotational motions in the camera frame. However, the distance to the camera remains the same over the sequence. We have also considered different lighting conditions to simulate day and night (where the background also changes) in order to prove the robustness of our methods whatever the situation. The full video result for the GPU-based approach is available at: https://youtu.be/H9XKDaQEVbU, and for the CPUbased approach at: https://youtu.be/Gr78pmRPKOA. For this first sequence, both methods have been able to perfectly track the target. The estimated trajectory, which is fairly the same for both methods, is accurate and smooth (see Figure 7). Furthermore, according to the collected ground truth data, our vision-based navigation provides precise pose estimation with an RMSE that is less than 1cm on the translation and about 2deg on the rotation (see Table 1). However, Table 1 also shows that the GPU-based solution is about 4 times more consuming than the CPU-based one (70ms per frame vs 290ms per frame), while considering twice less features (≈800 vs ≈420). This is mainly due to the fact that, for the GPU-based approach, we are using a 3D rendering engine with a much more complex 3D representation of the target. Finally, on Figures 8 and9, we can see the norm of the translational and rotational errors to the ground truth, which also shows where the algorithms encounter difficulties. Around frame 2200 for example, and for both approaches, the precision on the 3D pose estimation decreases a bit. The reason is that the panels of the target are orthogonally aligned with the focal axis of the camera (see first sequence, the satellite is here orbiting at 400km altitude, which is approximately the International Space Station (ISS) orbit altitude.

Let us note that this orbit corresponds to the one that will be used during the RemoveDEBRIS mission since the target will, in the first stages of the mission, be released from the ISS. Furthermore, during the mission, the CubeSat DS-2 will be ejected from the RemoveSat at very low velocity (≈2cm/s) and out of the orbit plan in order to optimize the relative trajectory for the vision-based navigation demonstration. In this experiment and as for the previous case, we considered different lightning conditions simulating day and night, which creates specularities on the satellite panels and consequently make more complex the navigation process. The CubeSat velocity has also been increased, with respect to the mission specifications in order to simulate larger inter-frame motions. Finally, this case also considers a non-fixed distance between the camera and the target so that the target appears sometimes very small in the image sequence with a very low visibility. Once again, this aims to demonstrate the robustness of the navigation process whatever the situation. The full video result for the GPU-based version is available at: https://youtu.be/AEJ6S7d9Y70, and at https://youtu.be/3nG4WQ9Tu6M for the CPU-based approach.

Regarding this second sequence, the CPU-based approach has not been able to correctly estimate the navigation data along the complete sequence (see Figure 12 and 13) while the GPU-based one has still been perfectly able to do it with an extremely good RMSE (0.0084 meters on the translation and radians on the rotation, see Table 2). By comparing Figures 13 and15 with Figures 10.c and 11.c around frame 1700, we are still facing the same issue as in the previous experiment. Indeed, for these frames, the panels of the target are orthogonally aligned with the focal axis of the camera, which creates incertitude on the orientation of the target. Furthermore, the target appears small with a low visibility (see Figure 10.c). As a result, the CPU-based approach has not been robust enough to correctly estimate the navigation data around frame 1700 (see Figure 13). Hopefully, as soon as the CubeSat has reappeared in the image with better conditions, this issue disappeared and the CPU-based approach has been able to recover correct results.

Note that in the presented approaches, no detection method has been implemented in parallel to the tracking. However, we could benefit from such algorithms to recover the target pose when tracking fails as, for example, in [START_REF] Petit | Robust visual detection and tracking of complex object: applications to space autonomous rendez-vous and proximity operations[END_REF][START_REF] Petit | 3D object pose detection using foregroundbackground segmentation[END_REF].

CONCLUSION

In this paper, we have presented two model-based tracking approaches that have been developed and tested on various experiments with different conditions. While the first method (which is available in the ViSP framework [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]) is only relying on a CPU, the second one is also taking benefit of using a low-end GPU. Furthermore, in addition to the edge-based and keypoint-based features also used in the first approach, the second one is also considering color-based features (that could be adapted to intensity-based features if needed). This has aimed to increase the robustness of the navigation process in order to face with complex situations, as shown for the second image sequence.

While considering space applications, the graphics processing unit might not be available and the central processing unit might not be powerful enough to run complex algorithms. However, when dealing with damaged or complex objects (which will probably be the case in the future of the RemoveDEBRIS mission), or when dealing with hard environment conditions as presented in the second experiment, the CPU version might not be robust enough and could eventually fail to estimate the navigation data. In those cases, the GPU-based approach remains a better choice.

Figure 1 :

 1 Figure 1: Model-based tracking principle considering the case of edge and keypoint features.

*

  edges that are hidden by other polygons. Neglecting the two above issues would result in adding extra outliers in the minimization process and consequently wrongly estimating the 3D pose of the camera.Once the visible edges of the model are projected in the image, each single line 𝑙 is sampled giving the set of initial image points { 𝐱 ' F } '34 M N . Each initial points is then tracked in order to determine the set of corresponding measured features { 𝐱 ' F * } '34 M N . As in[START_REF] Comport | Realtime markerless tracking for augmented reality: the virtual servoing framework[END_REF],

Figure 2 :

 2 Figure 2: Model-based tracking principle considering the case of color features.

  proposed methods have been validated on two synthetic image sequences where ground truth data are available to validate the quality of our visionbased navigation process. Both sequences simulate a CubeSat (DS-2) orbiting the Earth at different altitudes. The CubeSat DS-2 3D model has been generated by the Surrey Space Centre and provided by Airbus Defence and Space. The sequences have been generated via the 3D graphics and animation software Blender and the Earth has been textured using real satellite images available from NASA. Considered images dimensions are 640x480.

Figure 3 :Figure 4 :

 34 Figure 3: CubeSat DS-2 structure (image from[START_REF] Forshaw | An in-orbit active debris removal mission -REMOVEDEBRIS: Pre-Launch update[END_REF]).On Figure4.a, we can see the CAD model used for the GPU-based approach. It is a wireframe representation of the model generated by the Surrey Space Center. The model is fully detailed and composed by many triangles. On the other side, the CAD model that is used for the CPU-based approach has been manually re-designed from the detailed model in order to contain only visible edges in terms of gradient in the image. There is no triangle in this simplified version, and the visible polygons on Figure4.b directly represent the edges that will be used to define the edge-based and keypoint-based criterion errors (6) and[START_REF] Prisacariu | Pwp3d: Real-time segmentation and tracking of 3d objects[END_REF]. Models differences can also be seen on Figures 5.d and 6.d.

Figure 5 :Figure 6 :

 56 Figure 5: Vision-based navigation of a CubeSat DS-2for the first sequence of a geosynchronous orbit, using the GPU-based approach.

Figure 7 :

 7 Figure 7: 3D poses estimated for the first sequence of a geosynchronous orbit (translation in meters and rotation in radians).

Figures 5 .

 5 c & 6.c), which results in increasing the incertitude on the target orientation.Case 2: case of the ISS orbitIn the second sequence (see Figures10 and 11), we also use a CubeSat DS-2 target. With respect to the

Figure 8 :

 8 Figure 8 : Error to ground truth for the CPU-based approach on the first sequence.

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9 : Error to ground truth for the GPU-based approach on the first sequence.

Figure 12 :

 12 Figure 12: 3D poses estimated for the second sequence by using the CPU-based approach (translation in meters and rotation in radians).

Figure 13 :

 13 Figure 13: Error to ground truth for the CPU-based approach on the second sequence.

Figure 14 :

 14 Figure 14: 3D poses estimated for the second sequence by using the GPU-based approach (translation in meters and rotation in radians).

Figure 15 :

 15 Figure 15: Error to ground truth for the GPU-based approach on the second sequence.

Table 1 :

 1

First sequence: Root mean square errors (in meters for the translation and radians for the rotation), mean time and mean number of features per frame obtained for the estimated navigation values.

CPU-based MBT


Lesley-Ann A Duflot
email: lesley-ann.duflot@inria.fr-alexandre.krupa@inria.fr

Alexandre Krupa

Brahim Tamadazte
email: brahim.tamadazte@femto-st.fr-nicolas.andreff@femto-st.fr

Nicolas Andreff

Femto-St

Toward Ultrasound-based Visual Servoing using Shearlet Coefficients

This paper deals with the development of a new 6 degrees-of-freedom (DOF) vision-based controller for robotassisted medical applications. The main objective is to use visual information extracted from Ultrasounds (US) images to control an US probe held by the robot. Instead of the conventional use of visual features e.g., geometric features (points, lines, moments, etc.) in the visual control law design, the described method use the shearlet coefficients. More precisely, the time-variation of the coarsest level of shearlet decomposition coefficients are linked to the US probe (respectively to the robot) spatial velocity and then related to the task-function control law. The proposed control law was experimentally tested and validated using a realistic abdominal phantom. The obtained results demonstrated promising performances in terms of accuracy, repeatability, robustness and convergence behavior.

I. INTRODUCTION

2D and 3D US imaging devices, also known as echography systems, are widely used for non-invasive observations of inner organs i.e., diagnosis. Also, real-time US images can be used as visual feedback for guiding non/mini-invasive surgical procedures. Unfortunately, echography examination implies high demands on the surgeon since handeye-coordination becomes very mastery. Therefore, several robotic systems were developed to hold the US probe allowing thus remote examination of a patient. These systems are mainly based on master-slave robotic architectures allowing the clinician to tele-operate the probe motions [START_REF] Courreges | Clinical Trial and evaluation of a mobile, robotic tele-ultrasound system[END_REF]. Since few years ago, different solutions based on the visual servoing framework [START_REF] Chaumette | Visual servo control. part I: basic approaches[END_REF] were proposed to control the robot probe holder directly from information extracted from the US images. For example, in [START_REF] Krupa | Robotized teleechography: an assisting visibility tool to support expert diagnostic[END_REF] several autonomous control modes based on US image-based visual servoing were developed to provide different kinds of assistance during robotic US examination. The objectives of visual servoing during an echography examination can be summarized by: 1) maintaining visibility of an organ [START_REF] Abolmaesumi | Image-guided control of a robot for medical ultrasound[END_REF], 2) automatic positioning of the probe in such a way to retrieve a pre-operative image [START_REF] Mebarki | 2D ultrasound probe complete guidance by visual servoing using image moments[END_REF] and 3) automatic compensation of physiologic patient motion in such a way to actively stabilize the US image on a moving organ [START_REF] Nadeau | Intensity-based ultrasound visual servoing: modeling and validation with 2d and 3d probes[END_REF].

The efficiency of a visual servoing approach is strongly depended on the choice of the appropriate visual features. These features can be points [START_REF] Abolmaesumi | Image-guided control of a robot for medical ultrasound[END_REF], lines [START_REF] Hong | An ultrasound-driven needle-insertion robot for percutaneous cholecystostomy[END_REF], moments [START_REF] Mebarki | 2D ultrasound probe complete guidance by visual servoing using image moments[END_REF], etc. The success of these approaches depends highly on the ability to detect, match and track, in the US images, the visual features over time. In the literature, these techniques are referred to visual tracking. To overcome these limitations, original visual servoing approaches have been introduced which show that the design of an image-based visual controller can totally remove the visual tracking process [START_REF] Koichiro | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF], [START_REF] Tamadazte | A direct visual servoing scheme for automatic nanopositioning[END_REF]. Recently, in [START_REF] Nadeau | Intensity-based ultrasound visual servoing: modeling and validation with 2d and 3d probes[END_REF], the authors have proposed an alternative 6 DOF visual servoing scheme by taking as visual inputs directly the intensity of the US image pixels. This method has proven its reliability in terms of convergence domain, robustness and accuracy thanks to the redundancy of the information used as inputs to the controller.

This paper addresses another kind of featureless visual servoing undressed so far: the use of shearlet coefficients for control. Shearlet is a recent mathematics tool for signal (respectively, image) representation introduced in 2005 [START_REF] Easley | Sparse directional image representations using the discrete shearlet transform[END_REF] which can be considered as a natural extension of wavelet [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF]. It was proven that shearlet are more appropriate to represent anisotropic image features such as edges comparing to the wavelet. From this, several applications using shearlet representation have emerged: image inpainting [START_REF] King | Image inpainting: theoretical analysis and comparison of algorithms[END_REF], denoising and speckle filtering on US images [START_REF] Gupta | Enhancement of medical ultrasound images using multiscale discrete shearlet transform based thresholding[END_REF].

The main scientifc contribution of our work is the design of a visual servoing approach that use shearlet coefficients instead of pure photometric information to control the 6 DOF of a robotized US probe. One of the advantages using shearlet coefficients in the design of the control law is to obtain a set of noiseless and redundant visual signal for a more robust and accurate visual servoing. The latter was experimentally tested and validated using a 6 DOF US robotic system holder and a realist abdominal phantom. The validation tests were conducted under different conditions i.e., nominal and unfavorable. The obtained results were more than satisfactory especially in terms of convergence (smooth domain), robustness (robust to the unfavorable conditions), repeatability and accuracy (average error of 0.19mm and 0.23 • in translation and rotation axes, respectively).

In the remainder of this paper, Section II reviews the general basics of the shearlet representation of US images. The derived visual servoing controller is detailed in Section III while its experimental validation using a 6 DOF US robotic system holder and an abdominal phantom is discussed in Section IV.

II. SHEARLET TRANSFORM

The shearlet image representation used to design the proposed US-based visual controller constitutes an important mathematical tool that may allow development of promising new visual servoing approaches. Thereby, for a better understanding of the controller design methodology, it is obvious to introduce some basics of the shearlet transform, more precisely the cone-based shearlet system. As mentioned in [START_REF] Kutyniok | Introduction to shearlets[END_REF], the shearlet theory provides an efficient mathematical tool for image sparse representation including geometry and multiscale analysis. It can be considered as an extension of the wavelet transform, by increasing their directional sensitivity, to be more adapted for anisotropic image objects (e.g., edges, key-points, etc.).

Indeed, shearlet coefficients construction can be summarized by the reposition of three operators: dilatation D A α , shearing D S β and translation T t , applied on a generating waveform function φ (also called mother shearlet function). We invite the reader to refer to [START_REF] Kutyniok | Introduction to shearlets[END_REF], [START_REF] Kutyniok | Shearlab 3D: Faithful digital shearlet transforms based on compactly supported shearlets[END_REF] for more details.

Firstly, an operator D M is needed to compute the dilatation and the shearing function.

D M φ(x) = |det(M)| -1 2 φ(M -1 x) M ∈ GL d (R) (1)
with x ∈ R 2 , GL d (R) the group of d-dimensional invertible matrices, det(M) the determinant of the matrix M.

To obtain the dilatation function D A α we use this operator D M with the parameter M as a parabolic scaling matrix given by

A α = α 0 0 √ α α ∈ R + * (2)
Also, to obtain the shearing function D S β , we use the operator D M with the parameter M as

S β = 1 β 0 1 β ∈ R ( 3
)
where β is the slope that parametrizes the orientation of the generating function.

Finally, the translation function

T t ∈ L 2 (R 2 ) is applied using T t φ(x) = φ(x -t) (x,t) ∈ (R 2 ) 2 (4)
From this, the shearlet system associted to φ ∈ L 2 (R 2 ) is given by

SH(φ) = φ (α,β,t) = T t D A α D S β φ, α ∈ R, β ∈ R, t ∈ R 2
(5) By replacing (2), ( 3) and ( 4) in ( 5), we get

φ (α,β,t) (x) = α -3 4 φ A -1 α S -1 β (x -t) = α -3 4 φ 1 α -β α 0 1 √ 2 x -t (6)
Moreover, according to [START_REF] Kutyniok | Introduction to shearlets[END_REF] shearlet decomposition is possible if and only if φ satisfies the admissible condition:

R 2 | φ(ξ 1 , ξ 2 )| 2 ξ 2 1 dξ 1 dξ 2 < ∞ (ξ 1 , ξ 2 ) ∈ R 2 (7)
Thereby, the mother shearlet function φ ∈ L 2 (R 2 ) is defined such that :

φ(ξ) = φ(ξ 1 , ξ 2 ) = φ 1 (ξ 1 ) φ 2 ξ 2 ξ 1 (8)
where, φ is the Fourier transformation of φ, φ 1 is a wavelet function and φ 2 is a bump-like function. Then, the shearlet transform of the image i ∈ L 2 (R 2 ) is :

i → SH φ i(α, β,t) =< i, φ α,β,t > (9)
The effect of the parameters α and β on the frequency plane in the direction ξ 1 of the image i are shown in Figure 1, that will gives the horizontal coefficents. The same operations are made in the direction ξ 2 , which will gives the vertical coefficients. That's also important to notice that the low frequency box is isolated from the shearlet computation.

III. SHEARLET-BASED VISUAL SERVOING A. Basics of Visual Servoing

According to [START_REF] Chaumette | Visual servo control. part I: basic approaches[END_REF], the aim of a visual servoing is to control the motion of a robot in order to allow a set of visual signal s (s ∈ R k ) defining a robot pose r(t) ∈ SE(3) i.e., s = s r(t) to reach a set of desired signal s * (s * ∈ R k ) by minimizing a visual error given by e = s r(t) -s * (10)

The variation of s is linked to the velocity screw vector v = (v x v y v z ω x ω y ω z ) T of the visual sensor frame by . s = L s v where L s ∈ R 6×k is the interaction matrix. By taking into account [START_REF] Easley | Sparse directional image representations using the discrete shearlet transform[END_REF], we obtain the variation of the visual error e due to the visual sensor velocity such as . e = L s v -. s * . In order to obtain an exponential decrease of this visual error, the following classical control law is usually employed [

2] v = -λ L s + s(t) -s * ( 11
)
where λ is a positive gain and L s + is the Moore-Penrose pseudo-inverse of an estimation of the interaction matrix. In this study we choose to fix L s = L s * to be the interaction matrix L s * obtained at the desired position of the visual sensor as it is often done.

B. Control Law using Shearlet Coefficients

In this work, we propose to use the k shearlet coefficients of the coarsest level l = 5 (in the right of (Fig. 2 (b)) and in the bottom of (Fig. 2 (c)) as the set of visual features s for the image i, such as

s =< i, φ α,β,t > (l=5) α = 1, β = [-1, 0, 1] (12)
Since an analytical formulation of the interaction matrix related to the shearlet coefficients is difficult to derive, at this stage of development, we propose to numerically compute the constant interaction matrix L s * used in the control law [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF] thanks to an off-line initial robotic procedure that we performed at the desired location of the probe. The principle consists in applying successively and independently small relative displacements corresponding to 3 translations ∆x, ∆y, ∆z and 3 rotations ∆θ x , ∆θ y , ∆θ z along and around each axis of the Cartesian frame attached to the robot endeffector (as depicted in Fig. 7). The measured visual feature differences (scalar values) ∆s

C. Convergence domain

In order to judge the effectiveness of the proposed controller in term of convergence domain, we have computed the cost-function shape when varying the different robot DOF using the following relationship

C(s) = (s -s * ) (s -s * ) 1 2 (14)
As can be shown in Fig. 3, the cost-function C has an almost perfect convex shape where the global minimum is clearly identified.

IV. GROUND-TRUTH VALIDATION

To validate the proposed methods, we propose an automatic 6 DOF positioning task by controlling the 2D US robot holder motion. This scenario was performed in both simulation and experimental modes. The obtained results are presented and discussed below.

A. Simulation Results

Before the proposed approach can be validated on the experimental work-flow, it was tested using a developed C++ simulator which was very appropriate to validate the theoretical developments in different conditions of use without any risk of damaging the US probe. The simulator, miming a 2D US probe and an abdomen phantom, was implemented using the: Open-source Visual Servoing Platform (ViSP), Visualization ToolKit library (VTK) 1 and Shearlab library2 (available in MatLab version and implemented by ourselves in C++). An example of a US volume as well as a 2D US slice are shown in Fig. 4. This US volume is composed by 335 parallel 2D images of size 250 × 250 with a pixel size of 0.6mm × 0.6mm and interval of 0.3mm between each image (captured with a real 3D US Probe).

The test consists of performing a 6 DOF automatic positioning task using the shearlet-based visual controller. To do that, the positioning task consists of reaching a desired position r * ∈ SE(3) (i.e., desired image) from an arbitrary initial position r ∈ SE(3) (i.e., an initial image). At each US image acquisition, the shearlet coefficients are computed, more precisely the first 96 coarsest coefficients. The spatial velocities computed by the controller [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF] are applied to the virtual probe with a gain λ = 0.5 in order to reach

I di f f = (I -I * ) + 255 2 ( 15
)
The pose error obtained at the convergence using high resolution robot encoders supplied by the simulation software, was measured to be ∆e f inal (mm, deg) = 0.1, × , 0.2, 0.02, 0.18, 0.08 . Figure 6 shows the error decay for each DOF and the velocity evolution. After about thirty seconds, we can see that the system has converged.

B. Experimental Results [nominal conditions]

The proposed controller was tested and validated in experimental conditions using a realist testbench (Fig. 7) equipped with The objective of this test is to reproduce the simulation tests on a realist testbench closer to a clinical use of the US imaging system. For validation purpose, the US probe attached to the 6 DOF Viper robot is positioned to a reference location (preoperative image) where the desired shearlet components were measured and the related interaction matrix was numerically estimated, then without moving the phantom, we position randomly the probe at a new location which is considered as the initial position (intraoperative image). The preoperative image can be considered as a desired image acquired few days earlier and the intraoperative one captured during another examination. Thereby, the objective of our approach is to retrieve automatically and accurately the preoperative image during a second examination.

Furthermore, the y-translation corresponds to the contact direction between the US probe and the phantom. Indeed, it is advisable to be careful (low motion) along the y axis because of the risk of damage the probe and/or phantom. To tackle this problem, a force control is used to control this ytranslation instead of the visual controller (for more details, please refer to [START_REF] Nadeau | Intensity-based ultrasound visual servoing: modeling and validation with 2d and 3d probes[END_REF]) to maintain a desired contact force of 2N while the others 5 DOF are controlled by our shearlet-based visual servoing approach.

Figure 8 represents the case where the initial pose error between the initial US image I(r 0 ) and the desired US image I(r * ) was measured to be ∆e init (mm, deg) = (7, × , 8, 2,2, 2). The final positioning error (Fig. 9 (left), computed The positioning task was repeated several times using different initial positions (in order to avoid biased results). It was experimentally demonstrated that in each test the con-troller reaches successfully the desired position. TABLE I summarizes the initial and final errors obtained for each performed positioning task. Note that the y-translation is not taken into account since this axis was controlled by force control.

C. Experimental Results [unfavorable conditions]

In order to judge the robustness of the shearlet-based visual servoing, another experimental test was performed. During this test, the reference US image (Fig. 10(a)) was acquired in unfavorable conditions. Indeed, we learn the shearlet coefficients on the desired image blurred with a Gaussian filter of 21 pixels size. During the task, the current US images were captured in nominal conditions (Fig. 10(b)). The initial pose error is ∆e init (mm, deg) = (2, × , 1, 1, 5, 4). The positioning error decay in each DOF and evolution of the spatial velocity over the time are reported in Fig. 11(top left) and (top right). As can be underlined, the controller converges exponentially to the desired position. The obtained final error was measured to be ∆e f inal (mm, deg) = (0.219, × , 0.122, 0.14, 0.25, 0.25) very close to the ones obtained in the case of nominal conditions demonstrating the robustness of our method w.r.t. blurred images.

D. Towards Tracking using Shearlet-based Visual Servoing

The developed controller was also used to achieve tracking tasks. The latter consists in automatically moving the US probe to follow phantom motions. To do this, we use the same controller but set the desired image as being the initial one. Therefore the robotic system aims to automatically stabilize the US images through the compensation of the target external motions. The objective of this task is to compensate the physiological motions (e.g., breathing and heartbeat) during an US-based examination. For this experiment, we manually move the phantom along the 6 DOF with an amplitude of ∆A (mm, deg) = (20, 2, 2, 3, 17, 3).

As can be shown in (Fig. 12 and in the video accompanying this paper, our shearlet-based controller is able to efficiency compensate the applied external movements.

Fig. 1 .

 1 Fig. 1. Representation of a frequency plane tiled by the shearlet system. The horizontal cone is represented in blue.

Fig. 2 .

 2 Fig. 2. Shearlet coefficients (level l = 1...5) computed from the Lena photography (a). (b) the horizontal coefficients H SH(φ) and (c) the vertical ones V SH(φ).

Figure 2

 2 Figure2shows the representation of the shearlet coefficients of a 128 × 128 pixels image (e.g., Lena photography). These coefficients correspond to the different levels l = 1, . . . 5 (l = 1 corresponds to the highest resolution coefficients and l = 5 the lowest ones). Each level also contains a certain number of subimages, depending on the number of the chosen shearing directions.

  j ∆θ x , ∆ j s ∆θ y , ∆s j ∆θ z induced by these motions for each visual feature component j allows then to computed the j th row of the interaction matrix L s Finally, the whole interaction matrix is obtained by stacking all the rows L * j related to each visual feature component L s * = L * 1 , . . . , L * k and the control law (11) is applied with this numerically estimated constant interaction matrix L s * .

Fig. 3 .

 3 Fig. 3. (a) the shearlet-based visual servoing cost-function for x and y translation and (b) for the rotations around x, y and z axis.

Fig. 4 .

 4 Fig. 4. (a) 3D US volume provided by the developed simulator and (b) an extracted 2D US slice.

  Figure 5 depicts an image sequence acquired during the positioning task: Fig. 5(a)-(b) the initial and desired images, respectively, and Fig. 5(c) represents the initial image difference that is observed for an initial error of ∆e init (mm, deg) = (9, × , 5, 2, 2, 5) while Fig. 5(d) shows the image difference at the last recorded iteration, which shows the convergence of the control law.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Image sequence captured during the positioning task: (a) the initial image I(r 0 ), (b) the desired image I(r * ), (c) the difference I di f f between the initial and the desired image and (d) the difference between the final and the desired images.

• a 6

 6 DOF anthropomorphic robotic arm of type of Adept Viper s850; • a 2-5MHz 2D-US transducer (C60, Sonosite 180+) with a depth of 12cm providing 480 × 640 pixels US images with a pixel size of 0.35 × 0.36mm 2 ; • a 6 axes force sensor ATI Gamma SI-65-5; • a 3-GHz PC running under a Linux distribution in which where implemented the developed approaches in C++ and the control communication with the robot; • an abdominal phantom of type of AB-41900-030 Kyoto Kagaku-ABDFAN.

Fig. 7 .

 7 Fig. 7. Photography of the experimental set-up.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Image sequence captured during a positioning task: (a) the initial image I(r 0 ), (b) the desired image I(r * ), (c) the difference between the initial and the desired image and (d) the difference between the final and the desired images I di f f .

Fig. 10 .

 10 Fig. 10. Some snapshots captured during the positioning task [unfavorable conditions]: (a) the initial image I(r 0 ), (b) the desired image I(r * ), (c) the difference I di f f the initial image and the desired image and (d) the difference the final and the desired images.

Fig. 10

 10 Fig. 10 depicts some US images captured during the positioning task process. Fig. 10(a) and (b) show the initial image and the desired image (in poor quality) when Fig. 10(c) illustrates the difference between the initial and desired image I di f f and Fig. 10(d) the difference between the final and the desired image demonstrating the accurate convergence of the controller towards the desired position. For nominal or unfavorable conditions, the objective controller is to move the US probe (i.e., the holder robot) from the initial pose to the desired one by decaying smoothly the error in the shearlet domain toward zero Fig. 11 (bottom).The positioning error decay in each DOF and evolution of the spatial velocity over the time are reported in Fig.11(top left) and (top right). As can be underlined, the controller converges exponentially to the desired position. The obtained final error was measured to be ∆e f inal (mm, deg) = (0.219, × , 0.122, 0.14, 0.25, 0.25) very close to the ones obtained in the case of nominal conditions demonstrating the robustness of our method w.r.t. blurred images.

Fig. 11 .

 11 Fig. 11. Positioning task [top: experimental case, unfavorable conditions]: (left) error decay in each DOF, (right) velocity evolution in each DOF. [bottom: the evolution of the error norm e ]: (left) nominal conditions, (right) unfavorable conditions.

Fig. 12 .

 12 Fig. 12. Evolution of the error in the Shearlet domain for different tracking commandes.

TABLE I ROBUSTNESS

 I STUDY OF THE CONTROLLER ( T i (MM), R i (DEG), e 0 AND e f REPRESENT THE INITIAL AND FINAL ERROR, RESPECTIVELY.)

	e f and e 0	T x	T y	T z	R x	R y	R z
	e 0	1	×	2	5	5	5
	e f	0,41	×	0,09	0,06	0,13	0,14
	e 0	10	×	5	5	5	5
	e f	0,3	×	0,26	0,37	0,5	0,21
	e 0	7	×	8	2	2	2
	e f	0,16	×	0,09	0,09	0,22	0,07
	e 0	8	×	5	2	2	5
	e f	0,41	×	0,11	0,11	0,3	0,85
	e 0	4	×	5	5	5	5
	e f	0,3	×	0,06	0,1	0,22	0,08
	e 0	2	×	4	5	5	4
	e f	0,09	×	0,02	0,03	0,35	0,31
	µ(e)	0,28	×	0,10	0,13	0,29	0,28
	STD (e)	0,13	×	0,08	0,12	0,13	0,29

www.vtk.org

www.shearlab.org

V. CONCLUSION & FUTURE WORKS

In this paper was presented a new 6 DOF US imagebased visual servoing paradigm. Instead the conventional visual information in the visual control law design, our method makes use of the shearlet transform coefficients. This implied that the interaction matrix links directly the US images shearlet coarsest coefficients time-variation to the robot motion. The method was validated in both simulation and experimental using a 2D US probe attached to a 6 DOF anthropomorphic robotic arm. The obtained results in the cases of automatic probe positioning and phantom motion compensation demonstrated the efficiency of the proposed controller in terms of accuracy (average final positioning error of 0.19mm and 0.23 • in the translation and rotation, respectively), repeatability, robustness and convergence behavior.

Further work will be undertaken to improve the developed controller by providing an analytical formulation of the interaction matrix avoiding therefore the need of its off-line numerical estimation. We will also consider the problem of US images acquisition time and quality. This can be tackled by considering the Compressed Sensing techniques welladapted for medical images e.g., US, CT-scan, MRI, etc. The compressed sensing technique can also be taken into account during shearlet-based visual servoing.

VI. ACKNOWLEDGEMENTS

This work is partially supported by DESWEEP project funded by the Région de Bretagne and from the project NEMRO (ANR-14-CE17-0013-001) funded by the ANR, France. It is also performed in the framework of the Labex ACTION (ANR-11-LABEX-01-01).


MD Timothée Wallenhorst
email: timothee.wallenhorst@chu-rennes.fr

MD Mael Pagenault

MD Guillaume Bouguen

MD Laurent Siproudhis

MD, PhD Jean-François Bretagne

Small-bowel video capsule endoscopic findings of Cronkhite-Canada syndrome

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Cronkhite-Canada syndrome (CCS) is a rare, non-familial disorder characterized by multiple gastrointestinal polyps and ectodermal changes. This article presents the first small-bowel video sequences of CCS using video capsule endoscopy (VCE).

A 78-year-old female patient was admitted to our department with severe protein-losing enteropathy, alopecia, and onychodystrophy of the finger and toenails. The patient had no significant personal or family medical history. Endoscopy revealed multiple polyps of the gastric, duodenal, and colonic mucosa. Histological findings included cystic dilatation and gland distortions with inflammatory infiltration and stromal edema. Endoscopic and histological findings, as well as clinical symptoms, were typical of CCS. Further investigation by VCE (PillCam 2, Given Imaging, Israel) revealed a highly abnormal endoscopic appearance that has not been reported for other small-bowel diseases.

In the duodenum and jejunum, we observed whitish, excessively elongated and arborescent villi resembling dense underwater seagrass in the intestinal lumen. We also observed several herpes-like polyps, a feature that has been reported once. These abnormalities decreased as the capsule progressed along the small bowel, with reduction or absence of folds related to mucosal atrophy in the distal ileum.

In addition to gastroscopy and colonoscopy, small-bowel VCE may be useful for assessing the efficacy of CCS treatments.

Figure 1 :

 1 Figure 1: Jejunum video capsule endoscopy showing whitish, excessively elongated intestinal villi resembling dense underwater seagrass in the intestinal lumen.

This work did not receive any financial support.