Frédéric Cointault

L Journaux

Gilles Rabatel

C Germain

D Ooms

Marie-France F Destain

Nathalie Gorretta

G Grenier

O Lavialle

A Marin

Cointault Frédéric

Journaux Ludovic

Rabatel Gilles

Germain Christian

Ooms David

Grenier Gilbert

Lavialle Olivier

Marin Ambroise

Texture, Color and Frequential Proxy-Detection Image Processing for Crop Characterization in a Context of Precision Agriculture

published or not. The documents may come

Texture, color and frequential proxy-detection image processing for crop characterization in a context of precision agriculture

Introduction

The concept of precision agriculture consists to spatially manage crop management practices according to in-field variability. This concept is principally dedicated to variable-rate application of inputs such as nitrogen, seeds and phytosanitary products, allowing for a better yield management and reduction on the use of pesticides, herbicides … In this general context, the development of ICT techniques has allowed relevant progresses for Leaf Area Index (LAI) [START_REF] Richardson | Modeling Approaches to Estimate Effective Leaf Area Index from Aerial Discrete-Return LIDAR[END_REF], crop density [START_REF] Saeys | Estimation of the crop density of small grains using LiDAR sensors[END_REF], stress [START_REF] Zygielbaum | Non-destructive detection of water stress and estimation of relative water content in maize[END_REF] … Most of the tools used for Precision Farming utilizes optical and/or imaging sensors and dedicated treatments, in real time or not, and eventually combined to 3D plant growth modeling or disease development [START_REF] Fournier | ADEL-Wheat: A 3D architectural model of wheat development. Plant growth modeling, simulation, visualization and their applications[END_REF][START_REF] Robert | Coupling a 3D virtual wheat plant model with a Septoriatritici epidemic model: a new approach to investigate plantpathogen interactions linked to canopy architecture[END_REF]. To evaluate yields or to better define the appropriated periods for the spraying or fertilizer input, to detect crop, weeds, diseases …, the remote sensing imaging devices are often used to complete or replace embedded sensors onboard the agricultural machinery [START_REF] Aparicio | Spectral vegetation indices as non-destructive tools for determining durum wheat yield[END_REF]. Even if these tools provide sufficient accurate information, they get some drawbacks compared to "proxy-detection" optical sensors: resolution, easy-to-use tools, accessibility, cost, temporality, precision of the measurement … The use of specific image acquisition systems coupled to reliable image processing should allow for a reduction of working time, a lower work hardness and a reduction of the bias of the measurement according to the operator, or a better spatial sampling due to the rapidity of the image acquisition (instead of the use of remote sensing). The early evaluation of yield could allow farmers, for example, to adjust cultivation practices (e.g., last nitrogen (N) input), to organize harvest and storage logistics. The optimization of late N application could lead to significant improvements for the environment, one of the most important concerns that precision agriculture aims to address.

We propose in this chapter to explore the proxy-detection domain by focusing first on the development of robust image acquisition systems, and secondly on the use of image processing for different applications tied on one hand to wheat crop characterization, such as the detection and counting of wheat ears per m² (in a context of yield prediction) and the weed detection, and on the other hand to the evolution of seed development/germination performance of chicory achenes. Results of the different processing are presented in the last part just before a conclusion.

Image acquisition system

Image acquisition is the most important step for robust image processing. Indeed, the use of natural images involves some difficulties tied to outdoor conditions (lighting variations) and object complexity (the leaf area presents several planes, high contrast and lighting variations, scale variations, wheat growth stages ….). Photographic slides were taken at different growing stages of the wheatears.

Different solutions exist for image acquisition according to the vertical distance between the camera and the scene (1m, 2m or more), to the illumination used (natural or controlled light), to the kind of images (color, grey level, hyperspectral ...)

Figure 1 shows one solution based on specific boxes with opaque protection and controlled illumination by power-leds, for wheat ear counting, allowing to obtain color or grey level images. This device is constituted by a motorized rod installed on a tractor which allows a longitudinal displacement of a push-broom hyperspectral camera located 1 meter above the vegetation (figure 3). The images are acquired by a CCD monochrome camera (Hamamatsu C5405-70), with a resolution of 640 x 480 pixels and 256 gray levels. The system is enclosed within a black box to avoid interference from ambient light and is air-conditioned (Tectro TS27, PVG Int. B.V., Oss, The Netherlands) to maintain a temperature of 20°C. The lamp is optically isolated and white light did not escape inside the box. The blue filter transmitted an unexpected, small amount of infrared light between 770 and 900 nm (detected with spectrometer AVS-SD2000, Avantes, Eerbeek, The Netherlands).

Proxy-detection image processing for wheat ear counting

Proxy-detection imagery for yield prediction needs the determination of the three components of the yield: number of ears / m², number of grains per ear and the weight per thousand seeds. Wheat ear shapes are very variable, especially due to the different inclinations of wheatears. However, wheatears show a very rough surface whereas leaf surfaces are smoother (figure 4).

Leaf surface

Wheatear surface

Fig. 4. 3D surface of leaf and ear

Considering this specific roughness (figure 4), texture analysis can be used in order to discriminate wheatears, leaves and ground.

Texture feature extraction

In order to extract the texture features, a variety of methods have already been proposed in the literature and tested in practice. In this context, three families of texture feature extraction exist: statistical, spectral and structural methods. A successful texture classification or segmentation requires an efficient feature extraction methodology.

Two methods are discussed in this chapter. The first one is based on statistical analysis by cooccurrence matrices, and the second one is based on Fourier filtering.

Cooccurrence matrices

The principle of these matrices is to count how often pairs of grey levels occur in a digital image of texture. In fact, we take into account pairs of pixels that are located along a given direction defined by the angle  and separated by a fixed distance d.

Let

  The digital image : IJ f LL G   assigns a grey level to each pixel. Furthermore, we consider four directions ( = 0°, 45°, 90° and 135°) along which a pair of pixels can lie. Thereafter, for the distance d, we are able to obtain the unnormalized frequencies P( , ', ,d)  .

For instance, P( , ',0 , ) d    is defined as the cardinality of the set of pixel pairs having the following properties:

In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671

(, ) a n d (' , ' ) '

'= , =' fij fi j ii dj j          (1)
For the other angles, the corresponding frequencies are defined similarly. In this document, we consider the distance d = 1 (this means we consider neighboring pixels). Moreover, texture of the wheatears, and, more generally, of the pictures, are not oriented along a common direction. Thus, for each pair of grey levels we sum over the four frequencies (corresponding to the four angles) and we note: ' (,' , , 1 )

PP        (2) ' ' ' ,' P p P       
is the element of the normalized cooccurrence matrix.

Generally cooccurrence matrices are not used directly because of their size. The principle is to compute a set of measures from these matrices, such as those proposed by [START_REF] Haralick | Textural Features for Image Classification[END_REF].

  ' ' ,' log pp       (6)
The four coefficients were calculated for windows of 17x17 pixels surrounding each pixel. This window size is smaller than the ears' size and allows for reducing noise effects. The computational cost, the cooccurrence matrices were calculated over posterised pictures (32 grey levels), which avoids obtaining hollow matrices, considering the number of pixel in each window.

Table 1 shows the values of the four features for the 3 samples of the figure 5. The features computed from the cooccurrence matrices are used in a learning step in order to construct a discriminant function, called here textural function, which identifies efficiently the pixel belonging to the « ears » class.

Texture

Considering a trial set of windows, the use of a stepwise selection leads us to selects a subset of parameters (from the four ones describe above) to produce a good discrimination model.

At the first step, we compute the variance ratio F for each feature:

between group variance within group variance

F    =     2 .. . 2 . /( 1) /( ) ii i ik i i nxx g xx n g      (7)
g = 2 is the number of groups (one for ears, one for leaves and ground), n i is the number of windows for the group i

( i nn   ), . i
x is the mean of the group i, ..

x is the mean of the trial set, ik x is the value of the window k belonging to the group i.

The value of the variance ratio reflects the parameter's contribution to the discrimination if it is included in the textural function. The parameter that has the largest ratio is selected at first. Afterwards, at each step, the order of insertion is computed using the partial F-ratio as a measure of the importance of parameters not yet in the equation. As soon as the partial F-ratio related to the most recently entered parameter becomes non-significant the process is stopped.

In a second stage, a discriminant analysis with the selected parameters is computed. Because there are only two groups to discriminate, only one discriminant factor F D (textural function) is obtained.

In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671

Results from the cooccurrence matrices

In order to illustrate the efficiency of this textural function, the values of F D obtained on the image may be rescaled in the form of a gray level scale. For example, figure 6 shows an original image and the corresponding result after the first stage.

Results and discussion

In order to show its relevance, the Textural function is used in the energy function of a Maximum A Posteriori segmentation algorithm (Martinez de [START_REF] Martinez De Guerenu | Designing an automatic counting system for wheatears[END_REF]. Figure 7 shows separately the Wheatear and Non-wheatear regions resulting from this segmentation stage. Wheatear region pictures show false positive wheatear detecton (leaves or stems classified as wheatears). Non-wheatear region pictures show false negative detection (wheatears that have been forgotten).

Both kind of erroneous classification are circled on the pictures. For the whole set of sample images, 1 to 3 segmentation errors have been found for each image. As the images show an average of 20 wheatears, the error rate can be estimated around 10%. Most of the errors encountered are identified as badly lightened wheatears. This kind of error could be easily overcome by using artificial light instead of natural light conditions, or using color images instead of grey level images, as it will be shown in the next section.

Wheatear regions Non wheatear regions

Fourier filtering

To improve the previous works, one of the first solution has been to represent each image in a colour-texture hybrid space, and to use mathematical morphology tools (Cointault et al., 2008a) in order to extract and count the number of wheat ears.

Although this previous method gives satisfying results (around 10% of well detection, compared to manual counting), the hybrid space construction method is a supervised one and is limited by the objectivity of the operator. Moreover the statistical methods are dependent on the choice of the direction for the process, need an important computing time and are lightning-dependent.

The objective has been thus to propose new detection algorithms more rapids, robusts, and invariants according to image acquisition conditions, based on Fourier filtering and two dimensional discrete fast Fourier transform (FFT) [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF]. This approach includes three important steps: high-pass filtering, thresholding and cleaning of the image based on mathematical morphology operations [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF].

The three important steps

For the high pass Fourier filtering, a two dimensional FFT is performed on the target image (eq. 8). Based on the centered Fourier image, a high pass filter is applied in order to eliminate low frequencies in the FFT image (figure 8). The cut off frequency is empirically sized by a 10 pixels width disk mask as it is shown in the figure 8.

The thresholding of the resulting image is based on anInverse FFT and a predetermined threshold is applied in order to eliminate low pixel values which do not correspond to wheat objects (ground, leaf,…) (figure Fig. 9).

The cleaning step aims at eliminate remaining "non wheat" pixel groups, which are small and scattered. It lies on mathematical morphology operation and is performed with three sub steps:



First, a dilatation, which aims at making bigger and closer pixel groups in the image.



Then a blurring convolution with a Gaussian smoothing operator, followed by a thresholding, which eliminate too small groups of pixels. These small groups are considered as miss. This step makes smaller the pixel groups that correspond to wheat and then justify a third step.  Finally, another dilatation is performed which aims at regenerate size of pixel groups corresponding to wheat ears (figure 10).

In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671 In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671

Wheat ear counting estimation

The first image processing by thresholding and high pass Fourier filtering gives a binary image, composed of several pixel groups. Each group represents one or more wheat ears to be counted. In order to estimate this number we analyze the shape of each group considering two possible configurations:

 Pixel group presents a convex or nearly convex shape pattern: it is considered that only one ear is present in this kind of group (figure 11a).



Pixel group presents a concave shape pattern (figure 11b): it is consider that several ears are presents. In order to quantify the number of wheat ears, we estimate a shape index based on two features extracted from each pixel group:

1 c SS Q S    (9)
With S the surface and S C the convex hull surface of the pixel group (Figure 6c).

Coming from this index, a numberX of ears is attributed for the group, considering this interval:

(

1 ) / 1 0 / 1 0 XQ X   (10)
Another approach, based on skeleton analysis can be found in [START_REF] Germain | Non Destructive Counting of Wheatear with Picture Analysis[END_REF].

Results and discussion

A fast visual observation (figure 12) shows that only a small amount of groups corresponds to non-wheat things in image. More precisely, the biggest wheat ears are well detected and well separated from surrounded leaves.

Ears that lie the nearest of the ground, that are partially hidden or that are a little bit over exposed in the image are not well detected. Small amount of very big leaves also remain after cleaning step.

In order to test our image processing we performed algorithm on 40 images and compare the results with the mean of manual counting done by several experts. For the example, five images sample have been randomly chosen (table 2). Within the image set, a high variability of ear's number can be observed. With most of the images, high pass Fourier filtering method returns slightly higher counts. Based on the whole image sample, absolute difference between manual counting and image processing counting is contained under the value of 6%. The mean error obtained is 4,02% for this set of images.

The development of a high pass Fourier filtering approach method aims at creating an easily usable and adaptable method while obtaining a best wheat ear detection. Therefore, it is essential to compare this approach with previously used methods such as hybrid space method [START_REF] Cointault | In-Field Wheat ears Counting Using Color-Texture Image Analysis[END_REF]. For a visual comparison, high pass Fourier filtering has been applied onto images that have been previously treated with hybrid approach (figure 13).

In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671

Fourier Hybrid Fig. 13. Visual comparison between Fourier and Hybrid approach for two different images.

Visual comparisons with hybrid space results show that high pass Fourier filtering approach eliminate more non wheat objects within the images. Fourier approach separates more efficiently ear groups. Calculation time for Fourier approach is few sec per image while its few minutes for hybrid space, with the same operating system.

High pass Fourier filtering gives global satisfying results. Although a close range of settings has been determined, inverse FFT remains a parameter that has to be adjusted according to input image. An empirical value has been found and gives good results for most of images but it could be optimized with an automatic threshold selection such as k-means methods (MacQueen, 1967).

In the context of wheat detection, it has been observed that some ear objects are eliminated after cleaning step. These non-detections mainly correspond to near ground ears or ears massively hidden under leaves, hence, it should be relativized as too low ears may have development problems and may not be considerate within wheat yield. Ears that are located in over or under exposed part of image are not well detected but it is not due to the algorithm but to the quality of the acquisition, which is limited by the natural conditions. Small amount of very big leaves also remain after cleaning step and eliminate these artifacts constitute a further axis of development, including shape analysis in cleaning step.

In the context of wheat ear counting, it is observed that counting error percentage decrease with number of ears in images, hence, best results may be obtained with images representing more area and yet, more ears. Actually, worst error, 5,56%, is obtained with only 36 wheat ears in the image.

It is important to note that in most cases, Fourier approach returns slightly higher counts than manually counts. It should be due to missing detection, such as remaining leaves or over exposed part of images. Counts should be more precise with the including of shape analysis in cleaning step. In all cases, whatever the method used, the only way to obtain a right detection is to use 3-Dimensional information.

In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671

Proxy-detection image processing for weed-wheat crop discrimination by hyperspectral imagery

In the domain of weed control, since the environmental and economical stakes are particularly great, the herbicides are largely spread in order to assure sufficient yields for the whole field whatever the infestation rate.

The reason is essentially a technological one. Even if some low-cost devices are currently available to assure a localised spraying of herbicides on bare soil (vegetation detection by photoelectric cells), no commercial product allows a reliable and localised post-emergence treatment. Indeed, a such apparatus needs a sophisticated perception system, based on digital vision and allowing to distinguish between weeds and crops. The identification of varieties inside the vegetation is nowadays the principal lock to the development of localized weed control.

The corresponding research are numerous and can be divided into two main approaches [START_REF] Slaughter | Autonomous robotic weed control systems : A review[END_REF]:

 the spectral approach, in which the plant reflectance is the main parameter, using hyper-or multispectral images [START_REF] Feyaerts | Multi-spectral vision system for weed detection[END_REF][START_REF] Vrindts | Weed Detection Using Canopy Reflection[END_REF][START_REF] Vrindts | Weed Detection Using Canopy Reflection[END_REF]. The difficulty is then to propose spectral differences sufficiently robusts to lightning conditions  the spatial approach, based on spatial criteria such as plant morphology [START_REF] Chi | Leaf shape modeling and analysis using geometric descriptors derived from Bezier curves[END_REF][START_REF] Manh | Weed Leaf Image Segmentation by Deformable Templates[END_REF], plant texture [START_REF] Burks | Classification of weed species using color texture features and discriminant analysis[END_REF] … The main difficulty is tied here to the complexity and the natural variability of the scenes.

The study proposed in this chapter is tied to the previous approach: hyperspectral images of the wheat crop are acquired during the weed control period, and associated to specific preprocessing to avoid illumination conditions. After, the possible spectral discrimination between wheat and dicotyledonous weeds by means of chimiometric tools has been evaluated.

Spectral pre-processing

The images acquired by the camera (after a first internal processing taking into account the spectral sensitivity of the sensor) are images of luminance, which by definition depend on both the reflectance of the scene and lighting conditions. This is why a reference surface (gray ceramic), for which the reflectance Rc was measured in the laboratory, was systematically placed in each scene (figure 14). The average luminance Lc observed in the image to the reference can then be corrected reflectance for each pixel of vegetation :

Rf = Rc * (Lf/Lc) ( 11
)
where Lf is the observed luminance for that pixel.

It is important to note that the Rf value obtained is a "apparent" reflectance. Other phenomena are to be considered, namely:



The inclination of the observed leaf from the direction of the solar source, introducing a multiplicative factor of the light, and hence Rf  the possibility of specular reflection, adding a component to the spectrum colorless, and resulting in an additive constant of Rf.

In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671 These two phenomena are taken into account by applying a standard-centering operation (SNV, or "standard normal variate") on each spectrum of the image [START_REF] Vigneau | Potential of field hyperspectral imaging as a non destructive method to assess leaf nitrogen content in Wheat[END_REF].

Chemometric model

Because hyperspectral images include hundreds of contiguous spectral bands, they have a very high potential for spectral discrimination, compared to classical colour or multispectral images. Unfortunately, in a classification context, these kinds of high-dimensional datasets are difficult to handle and tend to suffer from the problem of the "curse of dimensionality", well known as "Hughes phenomenon" [START_REF] Hughes | On the mean accuracy of statistical pattern recognizers[END_REF]) which causes inaccurate classification (figure 14).

Fig. 14. Illustration of Hughes phenomenon

In this context, it is important to note that some news techniques already exist to overcome this problem. Among all possible methods, many recent works have focused on the area of dimensionnality reduction approaches in order to reduce this highdimensional datasets while maintaining the relevance of the information contained in the signal. [START_REF] Journaux | Texture Classification with Generalized Fourier Descriptors in Dimensionality In Agricultural Science[END_REF]. One of these approaches is the Partial Least Square Regression (PLS-R) , which builds a low-dimension subspace by determining a set of spectral "latent variables". Compared to other reduction methods such as Principal Component Analysis (PCA), the PLS takes into account both inputs and outputs to build its subspace, leading to better performances. The PLS-DA (PLS Discriminant Analysis) is an adaptation of PLS-R to discrimination problems.

In our case, a total of 335 spectra was increased step by step on a first reflectance image (figure 15) in three categories, namely wheat (157 spectra), dicotyledons (60 spectra) and ground (118 spectra).

A PLS-DA model was then determined on this data set by cross-validation using commercial software in chemometrics (The Unscrambler v9.7, CAMO Software AS, Oslo, Norway).

The resulting discriminant model, which involves 8 latent variables, was then exported and applied to all pixels in an image test, using a dedicated software developed in C + +.

In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671

Fig. 15. Calibration image and sample positions (the reference ceramic can be seen on the left)

Results

Figure 16 shows the image test and the classification results. We observe an excellent discrimination between wheat and weeds whatever the local conditions of illumination (shadow), with the exception of one type of weed (red circle) which was not present in the calibration image, and was therefore not included in the model. These results show the very high potential of hyperspectral image processing to extract pertinent information for agronomic applications.

Proxy-detection image processing for seed development and germination performance of chicory achenes by chlorophyll fluorescence

Several aspects of seed quality are tested routinely to minimize the risk of sowing seedlots that do not have the capacity to produce the desired crop. Amongst these tests, seed germination is important since it represents the percentage of pure seeds that have the potential to produce established seedlings in the field. The rate of germination corresponds to the reciprocal of the time needed for a given germination percentage to be reached [START_REF] Halmer | Germination rate[END_REF]. Accurate procedure of germination tests performed in laboratory is defined by the International Seed Testing Association (ISTA).

Besides these laboratory tests, seed processing lines includes cleaning machines that remove dust and waste material, and conditioning machines performing dimensional sizing, density sorting, and colour sorting. In conventional colour sorting, discoloured seeds are rejected, on basis of inspecting seeds individually to detect differences in reflected colour.

However, even after several sorting operations, some seed batches can contain a large proportion of viable seeds but still not sufficient for commercial use. These batches are lost because the viable and non-viable seeds cannot be separated using the conventional processing methods. The proportion of immature seeds in these lost batches is unknown. It would therefore be useful to provide a new, non-destructive method of distinguishing immature seeds from mature seeds in order to improve sorting processes.

In this context, the potential of fluorescence imaging [START_REF] Chen | Machine vision technology for agricultural applications[END_REF][START_REF] Nedbal | Chlorophyll fluorescence imaging of leaves and fruits[END_REF]) has been examined. The chlorophyll degrades during fruit ripening and the process of degradation was described by [START_REF] Barry | The stay-green revolution: Recent progress in deciphering the mechanisms of chlorophyll next term previous term degradation next term in higher plants[END_REF]. The chlorophyll is also a highly fluorescent molecule. Fluorescence occurs when some of the light absorbed by the chlorophylls is re-emitted at longer wavelength, typically between 650 and 750 nm. The fluorescent properties of chlorophylls have been used to evaluate the maturity of cabbage seeds (Brassica oleracea L.) [START_REF] Jalink | Chlorophyll fluorescence of Bassica oleracea seeds as a non-destructive marker for seed maturity and seed performance[END_REF][START_REF] Jalink | Seed chlorophyll content as an indicator for seed maturity and seed quality[END_REF]. The results showed that the magnitude of the chlorophyll fluorescence (CF) signal was inversely related to the quality of seeds. The relationship between the CF and germination performance was studied for tomato (Solanum lycopersicum L.) by [START_REF] Jalink | Seed chlorophyll content as an indicator for seed maturity and seed quality[END_REF] as in their previous study [START_REF] Jalink | Chlorophyll fluorescence of Bassica oleracea seeds as a non-destructive marker for seed maturity and seed performance[END_REF]. They concluded that seeds with an intermediate CF level were of the best quality, followed by seeds having a low CF signal. Seeds having a high CF signal were the worst. [START_REF] Konstantinova | Chlorophyll fluorescence sorting as a method for improvement of barley (Hordeum vulgare L.) seed health and germination[END_REF] measured the CF of barley grains (Hordeum vulgare L.) with a SeedScan I Laser Sorter (Satake, Stafford, TX. USA), using the principle developed by [START_REF] Jalink | Chlorophyll fluorescence of Bassica oleracea seeds as a non-destructive marker for seed maturity and seed performance[END_REF] but including a laser light source instead of a LED. They concluded that sorting a barley seed lot into six subsamples varying in CF values resulted in an optimal quality for the subsamples with low and intermediate CF signals. [START_REF] Suhartanto | Chlorophyll in tomato seeds: marker for seed performance?[END_REF] thoroughly described the relationships between the fruit CF, seeds CF and germination performance of tomato.

Image analysis

A specific image analysis code was developed with the GNU Octave language [START_REF] Ooms | Evaluation of chicory seeds maturity by chlorophyll fluorescence imaging[END_REF]. After applying background correction (the fluorescence values were divided by the reflectance signal of paper and multiplied by 60), the images were segmented and images of individual seeds were created, each of them being rotated along the main axis of the seed. The pappus side, which is brighter and larger than the radicle tip, was automatically detected on the basis of the mean width of the left half, its mean fluorescence intensity, the right half width and the right half fluorescence intensity. The accuracy of the detection was greater than 98%. The image was thereafter divided into the "pericarp zone" (Pe, 77% of the seed length) and the "pappus zone" (Pa, 23% of the length). The value of 77% was a compromise based on the observation of 100 seed images. The mean fluorescence values of the two zones were recorded for data analysis. The measurement system and image analysis are summarised in figure 17.

In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671

Results

The observed period corresponds to the phase of reserve deposition in the seed described in [START_REF] Bewley | Seeds: physiology of development and germination[END_REF]. Figure 18 shows the evolution of CF during this phase on the stalk, the weight parameters (dry weight DW and water content WC) and the germination performance (GP and GR). The dry weight increased, while the water content was still high at the end of the observed period (> 45 % on the stalk at 44 DAF, while the WC of stored seeds is about 6 %).  the CF decreased during the filling phase;  the end of the filling phase corresponded to the physiological maturity, where themaximal germination percentage and vigour is attained [START_REF] Black | The Encyclopedia of Seeds[END_REF].

On the other hand, the efficiency of CF features may be negatively affected by the large variability of individual measurements (random differences between individuals).

Current work aims at estimating the correlations between the CF features, the weight parameters and the germination variables in outdoors and greenhouse environments, and to assess the added value of CF features in comparison to weight, size and density features to distinguish between viable and non-viable seeds using sorting simulations.

Conclusion

To predict wheat yield or to determine wheat growth stages, remote sensing is not the only solution. Proxy-detection systems allow to acquire high resolution images to be treated by robust algorithms such as high pass Fourier filtering. For example, to predict

In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671 yield, the number of wheat ears has to be determined and is the result of a detection and a counting step. Detection has been done by developing specific image acquisition system and by implementing algorithms like colour-texture hybrid space or more classical image processing based on Fourier filtering, according to the frequential information (redundant information) included in the images. In this chapter, we presented the high pass Fourier filtering technique which gives satisfying and robust wheat ear detection with lower computing time. Moreover, we have compared the detection results with those obtained by representation in a hybrid space. Even if satisfying results are obtained for this qualitative experiment, some improvements should be done such as including an automatic threshold determination after the inverse FFT and an efficient shape analysis in order to obtain a finer wheat ear detection and better artefacts elimination.

This work has also revealed the possibility of extending the use of pattern recognition techniques and textural feature analysis to other applications, including the automatic determination of the wheat growth stage, aiming at creating a decision tool for farmers. At this stage, it seems that we should include an agronomic validation in order to propose more specific model to help farmers. In the same context one perspective could be to couple proxy-detection with satellite or aerial images. Moreover using multispectral data would probably improve even more the efficiency of such an approach.

The study of hyperspectral imaging demonstrates the potential to separate weed-culture.

A more comprehensive study must now be conducted to assess the robustness and the spectral and spatial minimum necessary resolutions in the context of an operational implementation.

Most fruits and seeds have measurable levels of chlorophyll, respectively in their pericarp or testa, and this chlorophyll degrades with time. In the case of chicory, the commercial seed is a fruit which cannot be hulled and the seed is not observable directly. The low amount of 10 chlorophyll implies the use of a highly sensitive device and because of the presence of a distinct pappus, the imaging of fluorescence is favourable. It would be profitable to identify the factors (variety, season, climate, hydric stress, etc.) influencing the evolution of CF to predict the characteristics of the decrease of CF and to predict if its decrease is always concomitant with the increase of the germination performance.

Finally, combination of several image acquisition systems should give more interesting results and open the door of the detection of other important crop characteristics.

Fig. 1 .

 1 Fig. 1. Example of image acquisition system (left) and its specific illumination by Led (right)Figure2shows some examples of color images of wheat crop at three different growth stages with controlled illumination.

Fig. 3 .

 3 Fig. 3. Hyperspectral acquisition device Other imaging system is based on chlorophyll fluorescence. A xenon light source (Hamamatsu Lightingcure L8222, model LC5, 150 W) is used to produce a white light which passes through an interference filter (03FIB002, Melles Griot, Carlsbad, USA) with a central wavelength 410 nm and width at half-maximum 80 nm to excite the chlorophyll-a. The light is conducted to the seeds by an optical fiber. The blue light issued from the filter is absorbed by the seeds chlorophylls, which resulted in fluorescence emission. A high-pass filter (665 nm, 03FCG107, Melles Griot 03FCG107) ensured the selection of the fluorescent signal from the blue light reflected by the object.

Figure 5

 5 Figure5shows three samples of ear ,soil and leaf textures: wheatearstexture reveals many transitions between very different grey tones. For leaves and soil the transitions are more progressive.

Fig. 6 .

 6 Fig. 6. Example of the textural functionF D (right) for a given wheatear image (left)

Fig. 7 .

 7 Fig. 7. Wheatear and Non-wheatear region results for MAP segmentation

Fig. 8 .Fig. 9 .

 89 Fig. 8. Example of wheat image (a), (b) its associated FFT projection, (c) cut off disk, (d) zoom of cut off disk.

  Fig. 11. Pixel groups with nearly convex (a) and concave shapes (b), convex hull of a concave shaped pixel group (c).

Fig. 16 .

 16 Fig. 16. Test image and classification results (light gray: ground, red: weed, green: wheat)

Fig. 17 .

 17 Fig. 17. The components of the chlorophyll imaging device and the main steps of image analysis.(1) Xenon lamp with blue optical bandpass filter (370 -450 nm), (2) optical fiber, (3) cone of blue light, (4) seeds from one capitulum, (5) blue light reflected by the seeds, (6) chlorophyll fluorescence (650 -730 nm), (7) highpass optical filter (665 nm), (8) CCD camera with zoom objective, (9) connection to the computer, (10) raw image, (11) background correction, then creation of individual images, (12) automated detection of the pappus using dedicated software and estimation of the levels of fluorescence FPER and FPAP[START_REF] Ooms | Evaluation of chicory seeds maturity by chlorophyll fluorescence imaging[END_REF].

Fig. 18 .

 18 Fig. 18. Evolution of chlorophyll fluorescence (FPER: pericarp, FPAP: pappus), dry weight, water content at harvest, germination percentage and germination rate of chicory seeds at each maturation duration on the stalk from 16 to 44 days after flowering. Means with confidence intervals of the means. The two following facts are in favour of the use of CF features for the differentiation of immature chicory seeds from mature ones, and as indicators of seed vigour:

Table 1

 1

		contrast	energy	entropy	inverse different moment
	ears	8.1	0.04	4.9	0.3
	leaves	2.5	0.22	3.6	0.6
	soil	0.3	0.94	2.0	0.8

. Values of the four Haralick's parameters associated with the three textural patterns

Table 2 .

 2 In Agricultural Science, InTech ed., 2012, 252 p., ISBN 9789535105671 Specifications of five experimental drying runs used for validation.

		Original images	Fourier results
	Image	Counting	Difference (%)
		Manual	image processing
	1	139	142	2,11
	2	36	34	-5,88
	3	90,5	94	3,72
	4	116,5	122	4,51
	5	136,5	142	3,87

Fig. 12. Results of wheat ear detection using high pass Fourier Filtering.


Switching controller for efficient IBVS

Guillaume Allibert, Estelle Courtial

Abstract-Image-based visual servoing (IBVS) is now recognized to be an efficient and robust control strategy to guide robots using only visual data. Classical IBVS is commonly based on the Cartesian coordinates of points in the image. Although the convergence of the visual features to the desired ones is generally achieved, classical IBVS can lead to unnecessary displacements of the camera, such as the camera retreat problem in the case of a pure rotation around the optical axis. In contrast, IBVS based on the polar coordinates of points, is well adapted to carry out rotations around the optical axis but less adapted to manage translations. To take advantage of the benefits of each approach, we propose a new strategy for visual servoing based on a switching controller. The controller switches between Cartesian-based and polar-based approaches thanks to a switching signal provided by a decision maker. With the proposed controller, unnecessary 3D displacements are minimized without any 3D reconstruction and visibility constraints can be taken into account. The local stability of the closed-loop switching system is proved. A comparison with classical controllers and advanced controllers is performed by simulations. This comparative study illustrates the effectiveness of the proposed controller in terms of displacement in the image space and in the 3D space.

I. INTRODUCTION

Visual servoing is a feedback strategy to control the movement of a robot by using only visual information [START_REF] Chaumette | Visual Servo Control, Part I: Basic Approaches[END_REF], [START_REF] Chaumette | Visual Servo Control, Part II: Advanced Approaches[END_REF]. Classically, two approaches are distinguished: image-based visual servoing (IBVS) and position-based visual servoing (PBVS). In IBVS, an error signal is measured in the image and the control objective consists in decreasing exponentially this signal to zero. Because it does not require 3D reconstruction, the IBVS strategy offers easy implementation and robustness in the face of measurement noise and model uncertainties. On the other hand, there is no control over the camera position. Unnecessary displacements that can lead to the system failure can be performed in the 3D space. Different IBVS approaches have been developed to avoid or reduce the useless displacement in the 3D space. The partitioned approach [START_REF] Corke | A new partitioned approach to imagebased visual servo control[END_REF] addressed this problem by decoupling and servoing the Z-axis rotational and translational components of the control. Two new image features (an angle and the square root of area) were computed. A potential function was also incorporated to keep the features in the image plane. In [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF], Chaumette proposed to use image moments to design a decoupled IBVS. When a combination of moments is well chosen, a satisfactory behavior of the features was obtained in the image plane as well as an adequate trajectory G. Allibert is with the laboratory I3S-CNRS/UNS, University of Nice Sophia Antipolis, France. Guillaume.Allibert@unice.fr E. Courtial is with the laboratory PRISME EA 4229, University of Orléans, France. Estelle.Courtial@univ-orleans.fr in the 3D space. To do this, the object and camera planes were parallel at the desired position and the object was a symmetrical object (but not a square object). The epipolar geometry was also used in hybrid methods [START_REF] Deguchi | Optimal motion control for image-based visual servoing by decoupling translation and rotation[END_REF], [START_REF] Malis | 2 1/2 D Visual Servoing[END_REF]. Based on the homography between the desired and initial images, these methods required coplanar feature points and an important computational load. An alternate approach [START_REF] Malis | Improving vision-based control using efficient second-order minimization techniques[END_REF] proposed different control schemes based on second-order minimization technique to avoid the camera advance/retreat problem. However, the estimation of the depth at the current instant was required. In PBVS, the camera position is estimated from the visual features extracted from the current image. If the pose is perfectly estimated, zeroing the error in the 3D-space provides a straight trajectory of the camera. There is no control in the image plane, which can result in the target getting out of the camera field of view. Several studies have addressed the complementarity of the two approaches [START_REF] Deng | Hybrid motion control and planning strategies for visual servoing[END_REF], [START_REF] Gans | Stable Visual Servoing Through Hybrid Switched-System Control[END_REF]. In [START_REF] Gans | Stable Visual Servoing Through Hybrid Switched-System Control[END_REF], a switched system using both IBVS and PBVS was developed. The switching was done when the image error or the pose error was considered too large. However, the approach remained susceptible to local minima and image singularities. In [START_REF] Deng | Hybrid motion control and planning strategies for visual servoing[END_REF], a hybrid control strategy consisting of a local switching control between IBVS and PBVS was developed in order to avoid image singularities, local minima and to take into account visibility constraints. While PBVS makes it possible to deal with local minima and image singularities, a 3D model of the target is required. Furthermore, PBVS is not robust to modeling and image processing errors when the features considered are points. In order to overcome the need of 3D information, the authors in [START_REF] Corke | Combining Cartesian and cylindrical coordinates in IBVS[END_REF] focused on the capabilities of IBVS by using two different parametrizations. The proposed hybrid visual servoing combined both Cartesian and polar image coordinates by stacking the respective interaction matrices. The induced trajectories were no longer straight lines in the Cartesian and polar image plane. The performance was a mean of performances, with no optimality. Furthermore visibility constraints could not be addressed. In this paper, we propose, not to combine but to switch between Cartesian and polar-based approaches in order to take advantage of the benefits of each approach. The most difficult problem is to determine when the switch should be done. A decision maker provides the switching signal according to the current value of the visual features and criteria. The latter are based on the prediction of the camera displacements up to the desired position. The switching controller, based on simple visual features like points (coplanar or not), minimizes unnecessary displacements of the camera without 3D reconstruction and can take into account visibility constraints.

The paper is organized as follows. Section II recalls the IBVS strategy respectively based on Cartesian coordinates, named IBV S C , and on polar coordinates, named IBV S P . The advantages and drawbacks of each controller are underlined. In Section III, the new controller is presented: the switching control architecture, the decision maker, the handling of visibility constraints and finally the stability analysis. In Section IV, numerous simulations show the efficiency of the proposed controller by comparing its performances to those obtained with classical controllers (IBV S C and IBV S P ) or with advanced controllers (controller using image moments [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF] or using second-order minimization techniques [START_REF] Malis | Improving vision-based control using efficient second-order minimization techniques[END_REF]). Conclusions are given in Section V.

II. IBVS IN CARTESIAN AND POLAR

COORDINATES The objective of IBVS is to regulate to zero the image error e between the current value of image features s and its desired value s ⋆ . The current value of image features is provided by measurements in the image plane which are usually the pixel coordinates of a set of n image points (n ≥ 3). In the following, theoretical developments are carried out for only one point.

A. In Cartesian coordinates

Consider P a 3-D point with coordinates (X,Y, Z) in the camera frame where Z is perpendicular to the image plane. Its projection in the image plane as a 2-D point s c makes it possible to write in the normalized case:

s c = (u, v) = (X/Z,Y /Z) (1)
where (u, v) are the Cartesian coordinates of an image point feature. The time derivative of (1) yields the fundamental relationship of IBVS:

ṡc = L s c τ (2)
where L s c ∈ R 2×6 is the interaction matrix and τ = [T x , T y , T z , w x , w y , w z ] ∈ R 6×1 is the spatial velocity of the camera. The interaction matrix based on Cartesian coordinates of one image point feature is given by [START_REF] Chaumette | Visual Servo Control, Part I: Basic Approaches[END_REF]:

L s c = -1 Z 0 u Z uv -(1 + u 2 ) v 0 -1 Z v Z 1 + v 2 -uv -u . (3)
Note that L s c depends on Z which is the depth of the point relative to the camera frame. This value must be estimated or approximated, and is noted Ẑ in the sequel. An IBVS controller is generally designed such that the image error exponentially decreases:

ėc = -λ e c with e c = s c -s ⋆ c and λ > 0, (4)
where s ⋆ c is the Cartesian coordinates of the desired image point. Combining ( 2) and ( 4), the resulting control law is then given by:

τ = -λ L s c + e c ( 5
)
where L s c + is the pseudo-inverse of the interaction matrix evaluated with the approximated depth Ẑ.

With this parametrization, the trajectories obtained in the image plane are straight lines when n = 3 and Z is perfectly known. Due to the particular expression of the interaction matrix (3), the IBV S C controller is well adapted to deal with translational displacements. Indeed, the camera translations along the X and Y axes are invariant to image point variations respectively along the u and v axes as can be seen in the first and second columns of (3). Nevertheless, in the case of pure rotation around the optical axis, the camera movement is composed of a rotational displacement but also of a translational displacement along the optical axis. This phenomenon is referred to as camera retreat [START_REF] Chaumette | Potential problems of stability and convergence in image-based and position-based visual servoing[END_REF]. This unnecessary movement is due to the coupling between the third and sixth columns in (3). To sum up, the IBV S C controller is generally efficient but in particular cases, it can generate unnecessary 3D displacements.

B. In polar coordinates

The coordinates of an image point feature can also be expressed in polar coordinates [START_REF] Chaumette | Visual Servo Control, Part II: Advanced Approaches[END_REF]:

s p = (ρ, θ ) = ( u 2 + v 2 , atan( v u )) = φ cp (s c ) ( 6
)
where ρ is the radius and θ is the angle. As previously, the following equations have ensued from the time derivative of ( 6):

ṡp = L s p τ, (7)
and the feedback control law:

τ = -λ L s p + e p . (8)
L s p ∈ R 2×6 is the interaction matrix based on polar coordinates. Its expression for one point is given by:

L s p = c Z s Z -ρ Z -(1 + ρ 2 )s (1 + ρ 2 )c 0 -s ρZ -c ρZ 0 -c ρ -s ρ 1 (9)
where c = cosθ and s = sinθ . L s p + is the pseudo-inverse of L s p evaluated with Ẑ. The image error is e p = s p ⊖ s ⋆ p where ⊖ is the modulo 2π subtraction for the angular component. Remark: a degeneracy occurs for the polar representation when a point is at the origin. In order to overcome this singularity, the coordinates (ρ, ρθ ) can be used [START_REF] Iwatsuki | A new formulation of visual servoing based on cylindrical coordinate system with shiftable origin[END_REF].

With this parametrization, the trajectories of the visual features are straight lines (n = 3 and the depth Z known) in the polar plane but obviously not in the Cartesian plane. As can be seen on the interaction matrix (9), the rotation around the optical axis Z is invariant to ρ (first row, last column) and directly proportional to θ (last row, last column). Furthermore, the translation T z is invariant to θ . This explains why the behavior obtained in the case of pure rotation around the optical axis is satisfactory. The camera movement is rotational only. On the other hand, translations along the X and Y axes require unnecessary displacements of the camera.

To take advantage of the benefits of each controller, a switching control strategy was developed.

III. SWITCHING CONTROLLER FOR IBVS

The switching control structure (Fig 1) makes it possible to commute from the IBV S C controller to IBV S P controller, and vice versa. The decision maker provides the switching signal according to output information s c and criteria.

A. Control law design

The aim of the proposed controller is to be able to switch between the two controllers in order to minimize or avoid unnecessary displacements of the camera. It goes without saying that if the movement required is a pure rotation around the optical axis, the controller chosen should be the IBV S P controller. Similarly, if the movement is a pure translation, the decision maker should choose the IBV S C controller. The problem is now to determine at which instant the controller should switch from one controller to the other for any movement. Since 2D errors in Cartesian and polar coordinates have not the same scale, they cannot be compared. The information common to both controllers can only be 3D information, namely the movement of the camera. The suitable controller should then guarantee that the desired position is reached from the current position while minimizing the movement of the camera. To determine the controller best suited to perform the displacement, a prediction of the camera displacement to be done until convergence is carried out in parallel with the two controllers. To perform the prediction of the camera displacement, a first step is the prediction of visual features [START_REF] Allibert | Predictive Control for Constrained Image-Based Visual Servoing[END_REF]. The prediction model is based on the dynamic equation (2) (resp. ( 7)) in a discrete-time version. With a classical first-order discretization scheme and considering Z to be approximated to Ẑ, the dynamic equations are written as:

s c (k + 1) = s c (k) + T e L s c (k) τ(k) ( 10
)
s p (k + 1) = s p (k) + T e L s p (k) τ(k) (11)
where T e is the sampling period.

Remark: Based on the expression of the interaction matrix, the prediction model has a local domain of validity. However, in practice, we encountered no problem of validity even for configurations that require large displacements.

Considering that one can roughly approximate the camera movement to the quantity T e T 2

x + T 2 y + T 2 z at each iteration, we can evaluate a view of the camera displacement required up to the desired position. The prediction procedure is given below for Cartesian coordinates. It is the same for polar coordinates: the subscript 'c' is changed to 'p'.

Prediction procedure:

1) At the current instant k, i = 0; 2) Measure s c (k + i) thanks to the image processing, 3) Update e c (k

+ i) = s c (k + i) -s ⋆ c and L s c (k+i) , 4) Compute τ(k + i) = -λ L s c + (k+i) e c (k + i), 5) Compute d c (k + i) = T e T 2 x(k+i) + T 2 y(k+i) + T 2 z(k+i) , 6) Predict s c (k + i + 1) = s c (k + i) + T e L s c (k+i) τ(k + i) 7) If e c (k + i + 1) = s c (k + i + 1) -s ⋆ c < β c then STOP else i = i + 1 and go to step 4.
β c is a small given constant and the corresponding time of convergence is noted (k + f c ) where f c is the number of sampling periods needed to achieve e c < β c . The prediction procedure is repeated at each iteration after new measurements. Therefore, we define the following two criteria:

J c = k+ f c ∑ j=k d c ( j), (12)
for the IBV S C controller and in the same way, for the IBV S P controller:

J p = k+ f p ∑ j=k d p ( j). ( 13
)
The decision maker should choose the IBV S P controller when J c > J p , and vice versa. Nevertheless, when the 2D errors are very small, the predicted 3D trajectories are very close for both controllers and consequently J c ∼ = J p . In order to avoid multiple switches of the controller, a hysteresis comparator has been implemented (see Fig. 2). If the IBV S C (resp. IBV S P ) controller was used at the previous instant, the IBV S P (resp. IBV S C ) controller will be used at the current time if and only if:

J p -J c > T -(resp. J p -J c < T + ) (14)
The parameter T is set by simulation. Finally, the decision maker which provides the switching signal is depicted on Figure 2. The prediction procedure and the decision making are repeated at each new sampling instant.

B. Visibility constraint

In the IBVS strategy, it is essential that the target remains in the camera field of view otherwise the feedback is no longer performed. This requirement is known as the visibility constraint. Due to the decoupling of the image errors, the IBV S C controller provides almost straight lines in the image plane (i.e. in the uv plane) if Z is estimated. The visibility of the target is then implicitly ensured. The problem of the target visibility can appear when the IBV S P controller is used. Indeed, it provides straight lines in the ρθ plane thus resulting in circular paths in the uv plane (Fig. 3). When the predicted visual features s p (k + 1) and/or s p (k + 2) are close to the safety zone defined in the image, the decision maker selects the IBV S C controller and keeps it during N iterations. Once the N iterations past, the standard procedure of prediction is restarted. As we will see later, N equal to 10 will be sufficient in practice.

C. Stability analysis

The stability of the closed-loop system with the classical feedback [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF] or [START_REF] Deguchi | Optimal motion control for image-based visual servoing by decoupling translation and rotation[END_REF] has been largely studied. The local asymptotic stability is proved in both cases [START_REF] Chaumette | Visual Servo Control, Part I: Basic Approaches[END_REF] under the assumption that the matrix L s L s + is positive definite (n = 3)

or positive semi-definite (n > 3. The matrix L s L s + is not full rank and local minima exist). Stability of the switching system can be difficult to prove especially under an arbitrary switching signal [START_REF] Liberzon | Switching in Systems and Control[END_REF]. This generally requires the stability of the controllers individually and also involves finding a Lyapunov function that is common to the two controllers.

Proposition 1: The closed-loop switching system (Fig. 1) is asymptotically stable for any arbitrary switching signal within a sufficiently small neighborhood of the origin.

Sketch of Proof: Firstly, a necessary condition is, above all, that each closed-loop subsystem with the IBV S C controller (resp. with the IBV S P controller) is individually stable. Then we have to find a candidate Lyapunov function common to both subsystems. If such a function exists then the closedloop system will be stable under any arbitrary switching signal. A mapping g can be defined between the image error in Cartesian coordinates and the image error in polar coordinates:

g : R 2 * n → R 2 * n e c → e p = g(e c ) = φ cp (s c ) -φ cp (s ⋆ c ), g(0) = 0. ( 15
)
Let us consider that V p is a common candidate function:

V p = 1 2 e T p e p . (16)
It is easy to show that V p = 0 if and only if e p = 0 or e c = 0 and V p = 1/2 e p 2 = 1/2 g(e c ) 2 > 0. V p is known to be a decreasing function for the IBV S P controller [START_REF] Chaumette | Visual Servo Control, Part I: Basic Approaches[END_REF]. We have to prove that it is also a decreasing function for the IBV S C controller. In the neighborhood of the origin e c = 0, the mapping g is injective and a Taylor expansion about e c = 0 gives:

g(e c ) = g(0) + dg de c | e c =0 e c + o 2 = Ge c + o 2 ( 17
)
where the jacobian G = dg/de c | e c =0 is full column rank and o 2 are terms of second order and higher. Combining ( 16) and ( 17), we get:

V p = 1 2 e T c G T Ge c + o 2 . ( 18
)
Considering ( 2) and ( 5), the time derivative of (18) gives:

Vp = -λ e T c G T GL s c L s c + e c + o 3 . ( 19
)
The function V p is a decreasing function in the neighborhood of the origin for the IBV S C controller if the matrix L s c L s c + is positive definite (or semi-definite). The candidate function V p is therefore a Lyapunov function common to both controllers. Consequently, the closed-loop switching system is stable under any arbitrary switching signal.

IV. SIMULATIONS

The control task consists in positioning a perspective freeflying camera with respect to a target composed of four coplanar points forming a rectangular object. The linear and angular velocities are limited to the following bounds: 25 cm/s for the translation speed and 0.25 rad/s for the rotation speed. The sampling period T e is 40 ms which corresponds to a usual camera (25 frames per second). To make the simulations more realistic, modeling errors of 20% are considered in the intrinsic camera parameters of the prediction models [START_REF] Gans | Stable Visual Servoing Through Hybrid Switched-System Control[END_REF] and [START_REF] Malis | 2 1/2 D Visual Servoing[END_REF]. Measurement noises are added to the outputs to simulate errors in the image processing.

A. Comparative study with classical controllers

In this section, the switching controller, noted IBV S S , is compared with respect to the IBV S C controller and the IBV S P controller. We consider that the depth is constant and equal to the depth at the desired position Z ⋆ . The following color code is used for the above simulations. For the image plane and the 3D plane, trajectories are plotted:

• in green for the IBV S C controller;

• in blue for the IBV S P controller;

• in red for the IBV S S controller. The switching signal is plotted in magenta on the control figure. When it is 0.3 (resp. -0.3), the controller IBV S P (resp. IBV S C ) is used.

Simulation 1: case of pure rotation. We first study a pure rotation (1.4 radian) around the optical axis of the camera. All errors in the image plane converge to zero (see Fig. 4(a)). The IBV S P controller is used throughout the simulation as expected (see Fig. 4(b)). Indeed, the IBV S C controller provides straight trajectories in the image but the corresponding 3D trajectory is made up of the rotation required around the Z axis and also of a translation along this axis. The latter is well-known to be useless as already mentioned in Section II. The trajectories in the image plane are circular in accordance with the 3D movement of the camera. The visual feature trajectories obtained with the switching controller are superposed to those obtained with the IBV S P controller. Simulation 2: case of pure translation. Translations along the X,Y, Z axes are tested. As already mentioned, the IBV S C controller provides an optimal trajectory of the camera in 3D space. The switching controller selects the controller IBV S C throughout the simulation (see Fig. 5(c)), thus proving that the controller chosen is the best controller in this case. The trajectories in the image plane and in the 3D space tend to straight lines from the initial to the desired pose (see Fig. 5(a) and 5(b)). No additional displacement is induced. Simulation 3: case of complex movement. The required movement of the camera is now composed of rotations and translations respectively around and along all the axes. We can see in Fig. 5(e) that the trajectories obtained with the IBV S S controller, in the 3D plane, are straightforward in comparison with those obtained with the IBV S C or the IBV S P controller. The camera displacement is clearly optimized: the useless displacements are minimized. The IBV S C is first selected for 0.5 sec and then the IBV S P is chosen, etc ... (see Fig. 5(f)). The trajectories obtained with the IBV S S controller in the image plane satisfy the best compromise between the convergence of the visual features and the straightforward displacement of the camera (see Fig. 5(d)).

To sum up, the decision maker provides a suitable switching signal whatever the camera motion to be done. The improvements with respect to classical controllers are significant: the 2D trajectory is still satisfactory while the 3D trajectory is improved by avoiding unnecessary displacement.

B. Comparative study with advanced controllers

In this section, the switching controller is compared with advanced controllers to show that its performances are similar to those of advanced controllers. Among the several approaches developed to design a decoupled IBVS as already mentioned in Section I, we have selected two of them: the method using second-order minimization techniques [START_REF] Malis | Improving vision-based control using efficient second-order minimization techniques[END_REF], noted IBV S PMJ , and the method using image moments [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF], noted IBV S IM .

Simulation 4: IBV S PMJ controller. In [START_REF] Malis | Improving vision-based control using efficient second-order minimization techniques[END_REF], two control laws were proposed. The simulation results presented here are performed by using the control law based on the pseudo-inverse of the mean of the Jacobians (PMJ). The two Jacobians are evaluated respectively at the desired position and at the current position, thus requiring the estimation of the depth Ẑ. In the case of pure rotation or translation, the simulations results of the two controllers are identical: due to the lack of place, simulation results are not shown. In the case of a displacement involving rotations and translations, the IBV S S controller provides the best result in both the 2D and 3D space (see Fig. 6(a) and 6(b)). Only one switch is needed (see Fig. 6(c)).

Simulation 5: IBV S IM controller. The IBV S IM controller is based on a combination of 6 image moments as described in [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF]. For a complex movement, we can see on Fig. 6(d) that the image trajectories of the IBV S IM controller are more direct. The behavior of IBV S S controller seems to satisfy a compromise between the trajectories in the 2D space and in the 3D space.

Although designed from simple visual features (points) and using conventional interaction matrices evaluated at the desired position, the effectiveness of the switching controller is similar to the advanced controllers's one.

C. Visibility constraint

The challenging difficulty of visual servoing is to keep the image features in the camera field of view. In order to illustrate the efficiency of the IBV S S controller, visibility constraints are taken into account in the relevant case of pure rotation (1.4 radian).

Simulation 6: Considering the image size to be ±0.3, a safety zone is delimited by bounds on u and v: u -= v -= -0.25 and u + = v + = 0.25. The convergence of the visual features to the desired position is still achieved while the visibility constraints are satisfied (see Fig. 7(a)). To do this, multiple switches were needed throughout the simulation (see Fig. 7(b)). While the IBV S P controller seems optimal from the 3D point of view for a pure rotation around the optical axis, the image trajectories can get out of the safety zone. A switch to the IBV S C controller ensures quasi-straight trajectories in the image and thus implicitly ensures the satisfaction of visibility constraints. Consequently, a slight translational displacement of the Although the control law generated by the controller IBV S S is not continuous, it is important to note that this control is applied as a reference input for low-level controllers servoing the joints of the robotic system. Therefore, the rapid changes are attenuated as was shown in practice in [START_REF] Gans | Stable Visual Servoing Through Hybrid Switched-System Control[END_REF]. In the worst cases, it is possible to filter the computed control before sending it to the low-level controllers. The average computational time devoted to the prediction of the camera movement for both controllers is about 20ms at each iteration and can be greatly reduced (with a C language) for a real time application.

V. CONCLUSIONS

In this paper, we have proposed a switching controller for IBVS tasks. The switching controller takes advantage of the benefits of Cartesian-based IBVS and polar-based IBVS. If the tools used (Cartesian IBVS, polar IBVS, switch, predictive approach) are well known, their combination is original and new. The switching approach is an alternate approach to design a decoupled IBVS in an easy way with point visual features and no 3D information. The visual features efficiently converge to the desired position while the 3D displacement of the camera is minimized. The decision maker provides the switching signal according to criteria based on the predicted displacements of the camera for the two controllers, without 3D reconstruction. The stability of the closed-loop switching system is locally proved for any arbitrary switching signal. The enhancement of the IBVS performances is illustrated through the simulation results. The robustness of the switching controller is checked in the presence of modeling errors and measurement noise. The handling of visibility constraints can be improved by introducing a potential function as in [START_REF] Corke | A new partitioned approach to imagebased visual servo control[END_REF]. The next step will be the implementation of the switching controller on a 6 ddl arm manipulator.

Fig. 1 .

 1 Fig. 1. Switching Control Structure

Fig. 2 .

 2 Fig. 2. Decision maker structure

Fig. 3 .

 3 Fig. 3. Visibility constraints

Fig. 4 .

 4 Fig. 4. Simulation 1: Case of a pure rotation

3 :

 3 Control inputs.

Fig. 5 .

 5 Fig. 5. Simulation 2: case of pure translation. Simulation 3: case of complex movement

  Simulation 5: 3D space.

  Simulation 5: Control inputs.

Fig. 6 .

 6 Fig. 6. Simulation 4: IBV S PMJ controller (in blue) vs IBV S S controller (in red). Simulation 5: IBV S IM controller (in green) vs IBV S S controller (in red).

Fig. 7 .

 7 Fig. 7. Simulation 6: Case of a pure rotation with visibility constraints


Hugo Talbot

Stéphanie Marchesseau

Christian Duriez

Hadrien Courtecuisse

Jatin Relan

Maxime Sermesant

Stéphane Cotin

Hervé Delingette

Interactive Electromechanical Model of the Heart for Patient-Specific Therapy Planning and Training using SOFA

Introduction

Simulating and personalizing electromechanical models of the heart can help in the training of cardiologists, as well as in planning the therapies and understanding the pathologies. We present here an interactive implementation of an organ-scale model of the heart with two applications: training for catheter ablation of cardiac arrhythmia and planning of cardiac resynchronisation therapy.

Cardiac arrhythmia is a very frequent pathology that comes from an abnormal electrical activity in the myocardium. Catheter thermo-ablation is a minimally invasive technique preventing such arrhythmias by removing the substrate or the trigger responsible for the cardiac dysfunction. Because of the ageing population, the number of patients suffering from arrhythmias is steadily increasing. Moreover this procedure requires highly skilled and experienced cardiologists. Therefore, the development of new training tools is required, and computerbased simulators can provide a good environment for this. To develop such a simulator, one important missing step is the ability to simulate the electrophysiology and the mechanical behaviour of the heart in real-time.

On the other hand, Cardiac Resynchronization Therapy (CRT) relieves Heart Failure (HF) symptoms with the implantation of a pacemaker that reduces heart dyssynchrony. However 30% of patients are non-responders to the therapy due the difficulties in configuring the pacemaker and selecting the patient. Many research groups (see for instance [START_REF] Sermesant | Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in CRT: A preliminary clinical validation[END_REF]) are investigating model-based solutions to tackle this significant issue for the clinical community. The hypothesis is that we can combine the anatomical and functional data to build patient-specific cardiac models that would have the potential to predict the therapy effects.

The contributions of this work are twofold. First, we developed an electrophysiological training simulator in SOFA which tackles the interactive issue in the context of cardiac arrhythmias. Coupled with this electrophysiology, we developed a mechanical model of the heart that can be personalized from MRI datasets.

Our simulations are based on the SOFA platform1 . SOFA is an open-source framework targeted at real-time simulation with an emphasis on medical simulation, mainly developed at Inria. A large choice of efficient solvers, hyperelastic or viscous material laws are already implemented in SOFA. Moreover, it enables interactivity during the simulation (pacing, surgery planning, ...) and gives a good trade-off between accuracy and computational efficiency.

Interactive Simulation of Cardiac Electrophysiology

We focus here on the training of cardiologists for the catheter-based ablation of cardiac tissues in the context of cardiac arrhythmias (see Fig. 1b). Our cardiac computations are based on the finite-element method. SOFA already enables the simulation of catheter navigation. In our simulation, we compute the navigation of a catheter inside the venous system from the femoral vein to the right atrium. In our work, we chose the Mitchell Schaeffer (MS) model [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF] to simulate electrophysiology since it provides a very good estimation of the transmembrane potential while computationally tractable for personalisation. We personalized the parameters of the MS model to fit patient-specific data [START_REF] Relan | Coupled personalization of cardiac electrophysiology models for prediction of ischaemic ventricular tachycardia[END_REF].

Targeting real-time computations is very challenging due to the stiffness of the equations involved. In order to use larger time steps, we first propose an adaptive parametrisation of the Mitchell-Schaeffer model to fit patient data. Second, a GPU implementation of the electrophysiology has been performed to decrease the computation time using recent methods [START_REF] Allard | Implicit fem solver on gpu for interactive deformation simulation[END_REF]. These optimizations allowed us to reach computations less than 3 times slower than real-time, thus making our simulation interactive. Same computations performed on CPU was at least 10 times slower than our GPU version.

Our simulation already enables the cardiology trainee to virtually stimulate any region of the heart. Such electrical stimulation are usually carried out to establish a definitive diagnosis before a thermo-ablation. Furthermore, the electrophysiological model can simulate thermo-ablation and the resulting electrophysiology. To do so, the factorized matrix of the system needs to be updated when ablating. The Bestel-Clément-Sorine [START_REF] Bestel | A biomechanical model of muscle contraction[END_REF] electromechanical model of the heart has been implemented in SOFA. It contains an hyperelastic component to account for the myocardium elasticity, implemented using optimized finite-element methods [START_REF] Marchesseau | Fast porous visco-hyperelastic soft tissue model for surgery simulation: application to liver surgery[END_REF], as well as an active part that deals with the contraction of the sarcomere. This active part uses the electrical potentials described previously as inputs for the mechanical contraction force. The model includes also viscous and elastic components to better characterize the fibre isotony during the deformation. Finally, the blood flow that circulates in the ventricles is also modelled globally using a constraint formulation. Efficient and interactive simulations were made possible thanks to the adoption of the simulation platform SOFA. Real-time simulations could not yet be reached but the computation time (about 10 minutes for one cardiac cycle) being reasonably small, a large set of personalization algorithms can be used on this implementation.

Patient-Specific Simulation of Cardiac Electromechanics

The authors developed a first calibration algorithm that enables to match the volume and pressure curves while estimating 7 parameters out of the 14 parameters of the complete Bestel-Clément-Sorine model. An example of the resulting deformation overlapped with the MR images during a cardiac cycle is shown Fig. 2.

Conclusion

In this paper, we presented our simulation environment based on the opensource framework SOFA for an interactive electromechanical model of the heart.

Our new performance results computing the cardiac electrophysiology have been detailed. Thanks to model optimizations and an efficient GPU implementation, we reached an interactive simulation close to real-time based on 3D patientspecific data. This first result is very encouraging regarding the development of our training simulator. Furthermore, first electromechanical personalization results were presented and proved that the now available model can fit MRI data in the context of patient-specfic therapy planning.

Such frameworks can open new possibilities in using computer models in the training of physicians as well as in the planning of therapies.

  (a) Simulated catheter navigation with synthetic X-Ray rendering. (b) Depolarization wave propagating inside the cardiac geometry.

Figure 1 :

 1 Figure 1: Simulation of catheter navigation and cardiac electrophysiology in SOFA.

Figure 2 :

 2 Figure 2: (Left) Short axis view of the simulated mesh overlaid on the MRI at different times of the cardiac cycle (heart period = 870 ms) on a pathological heart. (Right) Resulting Pressure-Volume loop.

http://www.sofa-framework.org