Benoit Kaufmann
email: benoit.kaufmann@gmail.com

Jean Louchet
email: jean.louchet@gmail.com

Evelyne Lutton
email: evelyne.lutton@inria.fr

Hand posture recognition using real-time artificial evolution

In this paper, we present a hand posture recognition system (configuration and position) we designed as part of a gestural man-machine interface. After a simple image preprocessing, the parameter space (corresponding to the configuration and spatial position of the user's hand) is directly explored using a population of points evolved via an Evolution Strategy. Giving the priority to exploring the parameter space rather than the image, is an alternative to the classical generalisation of the Hough Transform and allows to meet the real-time constraints of the project. The application is an Augmented Reality prototype for a long term exhibition at the Cité des Sciences, Paris. As it will be open to the general public, rather than using conventional peripherals like a mouse or a joystick, a more natural interface has been chosen, using a microcamera embedded into virtual reality goggles in order to exploit the images of the user's hand as input data and enable the user to manipulate virtual objects without any specific training.

Introduction

The work described in this paper is part of the REVES project 3 , whose aim is to create a prototype for a permanent exhibition at the Cité des Sciences et de l'Industrie de Paris 4 , called «Objectifs Terre : La Révolution des satellites»5 (objective: Earth). One of the components of this exhibition is devoted to artificial satellites and includes a virtual reality device that enables the user to interact with them.

The user's interface is based on the recognition of hand postures, using a fast generalised Hough transform operated using artificial evolution. The augmented reality device is described in Section 1.1 and the gestural interface in Section 1.2. We present the evolutionary Hough transform methology in Section 2, then describe two different implementations of hand model determination in Sections 3 et 4. Results are presented in Section 5 and the conclusions in Section 6.

The Augmented Reality device

The public will be staying around a large table. Above the table is a videoglobe 6 which displays an animated image of the Earth. Each user is provided with augmented reality "see through" goggles 7 (see figure 1) which allow to display the satellites that turn over the globe, or two-or three-dimensional information when the user is looking down at the table. To this end, the goggles are fitted with a miniature camera. Processing the images from this camera allows to position the synthetic images (e.g. the satellite images) displayed by the goggles relative to what the user can directly see through the goggles. It also allows the user to manipulate the objects displayed, by analysing the user's hand gestures.

Our contribution to this project was to design the algorithms that allow real-time recogniton of hand posture and position from the embedded camera and enable the user to interact with the virtual objects.

The googles.

Example of usage.

The gestural interface

In order to ensure natural interaction and allow the user to move whithin the available multimedia content, the interface should provide at least an object designation tool. The object may be e.g. a 3D object, an hypertext link, an icon or a menu entry. Selecting it allows to visualise richer information or move into different menu options. As this is assumed to replace the usual mouse click, we had to devise a simple communication language based on hand configuration. In addition to this, it is often useful to create a pointer which indicates where the system did locate the user's hand, in order to correct the eye-to-camera parallax and reduce the uncertainty that goes with the size of the detected patterns (e.g. pointed index or the whole hand). Then it is necessary to define two distinct gestures: one to move the pointer, and one to activate the object that has been pointed. The object is selected by changing the hand configuration. The advantage of the latter is the possibility to implement extra primitives such as rollover which allow access to other pieces of information by changing the pointer shape or the object colour. It is also possible to implement commands that move the display zone, or rotate a 3-D object. This could have been made possible through using the mouse buttons: a long pressure would switch from one shape to another one, and releasing the button would get back to the first shape. In fact we implemented this functionality by interpreting changes in the user's hand posture: closing and opening the hand is given the same semantics as pressing and releasing a mouse button. This improves the user's feeling to actually interact physically and grasp the object in order to turn or move it. Other display parameters (e.g. the scale) may be edited through creating cursors and moving them a similar way. It is also possible to relate the scale change to the apparent scale of the hand: a grasp plus a movement toward the user corresponds to zooming in, and reversely for zooming out.

The major problem to be solved here is the robust and real time detection of several hand gestures (at least, an open and a closed hand), without prior calibration or learning. To this end we developed an evolutionary Hough transform.

There exists an abundant literature about hand detection, but most of the proposed methods were not adapted to our conditions of use: for instance some methods are based on the detection of the entire body [START_REF] Darby | Activity classification for interactive game interfaces[END_REF][START_REF] Kim | Bimanual Hand Tracking[END_REF][START_REF] Carbini | Pointing gesture visual recognition by body feature detection and tracking[END_REF], or at least the user's foreharm [START_REF] Maimon | Hand Detection by Direct Convexity Estimation[END_REF] to deduce the position of the hands. This solution does not fit with our device because the camera is located on the user's face and therefore can only see the hand and not the rest of the body. Other methods only detect the hand but need a uniform [START_REF] Ionescu | Dynamic hand gesture recognition using the skeleton of the hand[END_REF] or at least a fixed background [START_REF] Zhao | An improved algorithm of hand gesture recognition under intricate background[END_REF]. Because the camera we are using is head-mounted, the user may turn their head, which makes the background change and show moving objects that could be wrongly detected as hands and disturb detection. For the same reason, motion-based [START_REF]Entertainment Computing -ICEC 2006, 5th International Conference[END_REF] and multi-camera [START_REF] Stødle | Gesture-based, touch-free multiuser gaming on wall-sized, high-resolution tiled displays[END_REF][START_REF] Kirishima | Real-time multiview recognition of human gestures by distributed image processing[END_REF] detection methods cannot be applied.

EvHough, evolutionary exploration of a parameter space

The Hough Transform [START_REF] Hough | Methods and means of recognizing complex patterns[END_REF] and its classical generalisations [START_REF] Maître | Un panorama de la transformation de Hough[END_REF] fundamentally scan the image looking for a certain type of local features. Each time such a feature has been found, it writes into the parameter space by reversing a model in order to increment the number of votes for the parameter values that are able to provide an explanation to the feature found. As the model is usually not injective, if n is the dimension of the set of parameters, each detected feature generates into the parameter space a variety with dimension n -1 : a curve if n = 2, a surface if n = 3, etc. This results into a very slow process with n = 3 and unrealistic processing times beyond.

The alternative proposed in [START_REF] Lutton | A Genetic Algorithm for the Detection of 2D Geometric Primitives in Images[END_REF] then in [START_REF] Louchet | From Hough to Darwin: an Individual Evolutionary Strategy applied to Artificial Vision[END_REF] consists in a direct, smart exploration of the parameter space, using a heuristics given by the Artificial Evolution paradigm. It does not involve any domain-specific knowledge (here, knowledge in image processing) except what has been expressed through the so-called "synthesis model" that allows to calculate the features corresponding to any point in the parameter space.

To this end, the algorithm creates a randomly initialised population of points in the parameter space, then evolves it using genetic operators: selection, mutation, crossover [START_REF] Rechenberg | Evolution Strategy : Nature's way of optimization[END_REF][START_REF] Baeck | A survey of evolution strategies[END_REF][START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF]. The selection process is based on a fitness function which, for each point in the parameter space, executes the model and calculates a value of similarity between the synthesised and the actual feature.

Thus, in the simplest case of sraight line detection, for each pair (ρ, θ), the evolutionary Hough algorithm will calculate the equation of the corresponding line, then evaluate e.g. the number of interest points on this line or the average contrast along it, in order to get the fitness value of the parameter pair (ρ, θ) in question.

The main interest of this is to open the possibility to work efficiently with high dimensional parameter spaces, where a direct Hough-style approach would fail. The main limitation is a classical one in the field of evolutionary algorithms: it remains uneasy to decide when the algorithm has to be terminated. However, in practice this is not a real issue: artificial evolution is naturally time-compliant and, unlike most classical image processing algorithms, may work on data that can be changed during processing. As shown below, this can be an interesting point when processing image sequences in real time.

3 Detection of hand models

Preliminary processing

The searched models are hand contours shapes, encoded as lists of points, in order to keep the number of points to be tested low. We thus apply a series of filters to the video sequences to be analysed.

Four convolution filters are first used to compute the gradient in 4 directions (vertical, horizontal and two diagonals). For each pixel of the input image, the largest of these 4 values then yields the gradient direction and value. The same information is also kept in memory for each point of the models.

In order to accelerate the computation and facilitate the optimisation process, a proximity function (equation 1) is computed at each frame, and gives for each x of the input image a measurement of its proximity to a contour [START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF].

prox(x) = max y∈image (grad y /dist x,y ) (1)
where x and y are two pixels of the input image, -grad y is the gradient value of pixel y, as defined earlier, -dist x,y is the distance between x and y according to equation 2 8 , with x i , for i = 1..8, the 8 immediate neighbours of x, as defined on figure 2. The gradient direction for each point of the input image is also updated in order to give the direction of the closest contour point.

dist x,x = 1 if x = y, dist x,y = min(dist x1,y + √ 2, dist x2,y + 1, dist x3,y + √ 2, (2)
dist x4,y + 1, dist x5,y + 1, dist x6,y + √ 2, dist x7,y + 1, dist x8,y + √ 2) x1 x2 x3 x4 x x5 x6 x7 x8
In the same time, the colour encoding of the input image (RGB) is transformed to Lab coordinates. In this colorimetric space, a skin colour model has been defined as a cylinder parallel to the luminance axis whose basis is circular in the space of chrominance (corresponding to an average range of chrominance coordinates of skins). This allows then to compute a low resolution binary image IndC, which reads which pixel (or corresponding square area of the input image) may correspond to a skin colour or not. Several cylinders can be defined according to different skin colours, in which case the IndC binary image will indicate if the corresponding pixel area corresponds to at least one of the colour cylinders.

Genome and similarity criterion

The search space, i.e. the space of all possible solutions, is a 5 dimension space, which corresponds to the 5 following parameters:

model for the hand model (see figure 4) -tx for the horizontal translation of the model, -ty for the vertical translation, -scale for the apparent scale of the model, which varies with the size of the user's hand and its distance to the camera, -rotate for the rotation of the hand with respect to the optical axis of the camera.

The other rotations are ignored, as we consider that the user designates objects with arms stretched, and hands quasi othogonal to the optical axis. This allow to restrict the search to 2D models of hands.

The similarity criterion (or f itness) is a combination of three independent criteria, f 1 , f 2 and f 3 .

f 1 gives a measure of the distance between the image and model contours, it is based on the use of the proximity image (see figure 3):

f 1 (S) = p∈contour p × prox(proj(p)) card(contour) (3)
S = [model, tx, ty, scale, rotate] is a point of the search space, proj(p) is the projection of p via the transformation defined by S and card(contour) is the number of contour pixels. Dividing by card(contour) allows the criterion to be independent of the number of contour pixels.

f 2 measures the correpondance of gradient directions between image area and model:

f 2 (S) = card(C) card(contour) (4)
card(C) is the number of contour pixels x whose direction corresponds to the model's direction, and card(contour) is the number of pixels of the contour. -f 3 uses the binary image IndC to verify if the interior of models contours has the right colour (this calculation is usually done on a lower scale image):

f 3 (S) = card(D) card(E) (5)
card(D) is the number of interior pixels of the model whose projection on the image corresponds to "true" on IndC and card(E) is the number of interior pixels of the model. Once again dividing by card(E) allows to have a measurement independent to the size of the model.

The resemblance criterion for a solution S is then:

f itness(S) = α 1 × f 1 (S) + α 2 × f 2 (s) + α 3 × f 3 (S) α 1 + α 2 + α 3 (6)
In our implementation, we choose α 1 = 10, α 2 = 10 and α 3 = 1.

Genetic engine

A steady state approach has been chosen for the EA engine, which means that each new individual immediately replaces an old unadapted one within the population: no generation synchronism is used. The individuals to be replaced are chosen using the tournament technique. Three genetic operators are used:

immigration: a new individual is randomly generated, mutation: a new individual is a small perturbation of an existing one, crossover: a new individual is created as a random weighted average of two existing individuals.

After an immigration or a mutation is done, two individuals are chosen in the population and the one who gets the lower fitness is replaced with the new individual (tournament of size 2). After a crossover, a tournament of size 3 is used. Steady state engine and tournament selection have been chosen in order to better maintain the diversity and limit the risk of premature convergence. This point is actually critical here, as we operate a small sized population in order to obtain real time performance. A generational process would necessitate fitness ranking which is time consuming. For the same reason, a fixed size population is used, in order to minimise time consuming memory allocations.

Additionnally, as this EA runs on an varying environment, the fitness function is recomputed at each new frame. We actually get an asynchronous algorithm, that uses video information at the time it is available. Finally, to compare with a generational EA, we measure the evolution in terms of "pseudo-generations", which corresponds here to the evaluation of a number of individuals equivalent to a generation gap (i.e. 40% of the population size, see section 5).

Experiments based on the previous algorithm yield good results, however in some conditions, some model representatives were lost in the population. As a consequence, the detection of a model change was unstable and slower. We thus preferred an encoding that naturally balances the number of representatives within the population, by not explicitly using the model any more: the fitness is then computed for each searched model. Each individual thus has now n fitness values, one for each of the n searched models. The resulting fitness is the maximal value of the n fitness:

f itness(i) = max(f itness 1 (i), ..., f itness n (i)) (7)
f itness(i) is the fitness value of individual i and f itness k (i) is the fitness value of individual i according to model k (k ∈ [1..n]), computed using equation 6.

For each frame, the detected model is the one which corresponds to the highest fitness, as soon as it is greater than a fixed threshold (0.3 for the results presented in the next section). To avoid unstabilities due to the definitions of a strict threshold, the detection decision is based on a double threshold and takes into account the model recognised at the previous frame:

if the largest fitness of the last population of the current frame is larger than threshold T , the corresponding model is detected, if the largest fitness is in [0.7 × T, T ], and corresponds to the detected model of the previous frame, it is maintained as "detected," else, nothing is detected.

Results

The results presented below have been obtained on a netBook MSI U100, equipped with a single core 1,60 GHz Intel Atom N270 processor, with 1 GB of memory, and using GNU/Linux Mandriva 2009. The images are captured with the integrated webcam (resolution of 160 × 120 pixels).

Examples of detection on a video sequence are displayed on figure 6, and the curves of figure 5 give a representation of the distribution of representatives of each of the searched models within the population. Each graph (one per searched model, see figure 4) presents three curves which are the maximal and average of fitness values, and the proportion of individuals whose fitness is at least equal to 90% of the maximal fitness. The horizontal line represents the model detection threshold (0.3). Vertical lines mark the transition between two successive frames. For each frame, its number, and the number of generations that were run on it, are displayed.

It can be noticed that the curve corresponding to the fitness max of model 1 (which is not present in this part of the video) never crosses the 0.3 threshold. On the contrary, one may notice on graphs of models 2 and 3, an alternance of detection of these models, and for each detection phase, a concentration of the population around the maximal values (green curve).

Other examples of detection are given on figure 7, in the experimental conditions described in section 1.1. A video sequence is also available on http://apis.saclay.inria.fr/twiki/bin/view/Apis/HandGestureRecognition

Conclusion

We showed that a real-time evolutionary algorithm can be run on a general-public device, using low computation power. The specific components of this evolution algorithm are a small population, an asynchronous steady-state genetic engine, that produces a variable number of generations between two frames, depending on the computation load of the whole system. An implementation on another computer would result in a different average number of generations between two frames. This characteristic yields an efficient capability to the algorithm, as it is able to exploit the available data as soon as they appear. Of course the number of generations between two frames condition the quality of the detection, and its reactivity to sudden changes. In extreme conditions, if the algorithm cannot obtain enough computation power, the population does not have enough time to converge between two frames. From the user's side, this results in late detections, and unability to follow rapid moves. The genetic encoding finally adopted (section 4) has been preferred for robustness and reactivity reasons. In case the genome carries the model, a delay is necessary in order to let representatives of a new model grow, while in the second version, this change can be immediate. Of course this solution is only efficient if the number of searched models is low, as fitness calculation is longer (n fitness computations for n models).

Additionnally, for the application described in section 1.1, a Kalman filter on the detected parameters for each frame allows to obtain a smoother detection. Finally, the display of a "mouse pointer" in the virtual environment (yellow mark in figure 7) gives functionalities similar to classical computer mouses.

Fig. 1 .

 1 Fig. 1. The interactive device (courtesy Zile Liu, Laster Technologies).

Fig. 2 .Fig. 3 .

 23 Fig. 2. Neighbourhood of x

Fig. 4 .

 4 Fig. 4. Hand models: the white pixels are contours, blue pixels are pixels where the colour model is evaluated, and black pixels are ignored.

Fig. 5 .

 5 Fig. 5. Variation of the population models distribution and fitness values, for a sequence where the 3 models of figure 4 are searched.

Fig. 6 .

 6 Fig. 6. Example of detection on a video sequence in natural light environment. Blue shapes show the detected model, red points correspond to the position (origin of the models) (x, y) of all individuals of the current population.

Fig. 7 .

 7 Fig. 7. Example of a detection on a video sequence in controlled light environment. A yellow marker represents the mouse position.

REVES ANR (French National Research Agency) contract No 2160(2007)(2008)(2009)

http://www.cite-sciences.fr/english

http://www.cite-sciences.fr/objectifs-terre

Developed by Agenium (http://www.agenium.eu/sections.php?op=listarticles&secid=40)

Developed by Laster Technologies (http://www.laster.fr/?lg=en)

This distance is an approximation of the euclidean distance, that is computed faster.

Acknowledgments:

The authors thank Sergio Alejandro Mota-Gutierrez, of Guanajuato University, Mexico, student of the Electrical Engineering Master Program, for his contribution to the colour based fitness computation.


Franck P Vidal

Jean Louchet

Jean-Marie Rocchisani

Évelyne Lutton

New genetic operators in the Fly algorithm: application to medical PET image reconstruction

This paper presents an evolutionary approach for image reconstruction in positron emission tomography (PET). Our reconstruction method is based on a cooperative coevolution strategy (also called Parisian evolution): the "fly algorithm". Each fly is a 3D point that mimics a positron emitter. The flies' position is progressively optimised using evolutionary computing to closely match the data measured by the imaging system. The performance of each fly is assessed using a "marginal evaluation" based on the positive or negative contribution of this fly to the performance of the population. Using this property, we propose a "thresholded-selection" method to replace the classical tournament method. A mitosis operator is also proposed. It is triggered to automatically increase the population size when the number of flies with negative fitness becomes too low.

Introduction

Image reconstruction in tomography is an inverse problem that is ill-posed: a solution does not necessarily exist (e.g. in extreme cases of excessive noise), and the solution may not be unique. This problem can be solved as an optimisation problem, and on such cases, evolutionary algorithms have been proven efficient in general, and in particular in medical imaging [START_REF] Bosman | Evolutionary algorithms for medical simulations: a case study in minimally-invasive vascular interventions[END_REF][START_REF] Cagnoni | Genetic algorithm-based interactive segmentation of 3D medical images[END_REF][START_REF] Völk | Multiple network CGP for the classification of mammograms[END_REF]. We focus here on positron emission tomography (PET) reconstruction in nuclear medicine.

Nuclear medicine appeared in the 1950's [START_REF] Badawi | Nuclear medicine[END_REF]. Its principle is to diagnose or treat a disease by administering to patients a radioactive substance (also called tracer) that is absorbed by tissue in proportion to some physiological process. When a pathology occurs, the metabolism most of the times increases: there are more molecules in the pathology area, i.e. the radioactivity also increases.

It is possible to reconstruct slices through the human body using methods similar to those used in conventional X-ray computed tomography [START_REF] Kak | Principles of computerized tomographic imaging[END_REF]. In nuclear medicine, this method makes use of a gamma emitter as radio-tracer. It is called Single-Photon Emission Computed Tomography (SPECT). The reconstruction allows to recover the 3D distribution of the tracer through the body.

The other main tomographic technique in nuclear medicine is PET. Here a positron emitter is used as radionuclide for labelling, rather than a single gamma emitter. Positrons are emitted with high energy (1 MeV). After interactions, a positron combines with an electron to form a positronium. Then the electron and positron pair is converted into radiations. It is the annihilation reaction, which generally produces two photons of 511 keV emitted in opposite directions. Taking advantages of this property, this radiation is detected in coincidence, i.e. using the difference in arrival times of the detected photons of each pair, and considering that each annihilation produces two photons emitted in exactly opposite directions. The line between the detectors that have been activated for a given pair of photons is called "line of response" (LOR). Prior to the reconstruction, the LOR data is often rebinned into a sinogram [START_REF] Fahey | Data acquisition in PET imaging[END_REF][START_REF] Lewitt | Overview of methods for image reconstruction from projections in emission computed tomography[END_REF]. This intermediate data representation corresponds to projection data that can be used by conventional tomographic reconstruction codes. A broad overview of reconstruction methods using projection data in nuclear medicine can be found in [START_REF] Lewitt | Overview of methods for image reconstruction from projections in emission computed tomography[END_REF][START_REF]Quantitative Analysis in Nuclear Medicine Imaging[END_REF].

The PET reconstruction methods are often divided into two classes: i) analytical methods, and ii) iterative statistical methods. Analytical methods are based on a continuous modelling and the reconstruction process consists of the inversion of measurement equations. The most frequently used is the filtered back-projection algorithm (FBP) [START_REF] Kak | Principles of computerized tomographic imaging[END_REF]. Statistical methods are based on iterative correction algorithms. These include the most widely used techniques in SPECT and PET, such as the maximum-likelihood expectation-maximization method (ML-EM) [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF] and its derivative, the ordered subset expectation-maximization algorithm (OS-EM) [START_REF] Hudson | Accelerated image reconstruction using ordered subsets of projection data[END_REF].

In a previous paper, we showed that a cooperative coevolution strategy (also called Parisian evolution) called "fly algorithm" [START_REF] Louchet | Stereo analysis using individual evolution strategy[END_REF] could be used in SPECT reconstruction [START_REF] Bousquet | Fully three-dimensional tomographic evolutionary reconstruction in nuclear medicine[END_REF]. Here, each fly corresponds to a 3D point that is emitting photons. The evolutionary algorithm is used to optimise the position of flies.

However, PET has taken over SPECT in routine clinical practice. Effort has therefore been made to propose an efficient evolutionary scheme that takes into account PET data acquisition principles [START_REF] Vidal | Artificial evolution for 3D PET reconstruction[END_REF][START_REF] Vidal | PET reconstruction using a cooperative coevolution strategy in LOR space[END_REF], but these were still restricted to low resolution PET scanners in 2D-mode. This paper describes our current research activities aimed at providing an effective method in both 2D or fully-3D mode, and it describes recent developments, such as i) the introduction of "thresholded-selection" replacing the traditional "tournament selection" and ii) taking advantage of the thresholded-selection to increase the population size (when the number of flies, whose fitness is negative, is too low), i.e. improve the statistics of the final image. The following section gives an overview of the methodology. The results and performance of our method using numerical phantoms are presented in Section 3. The paper ends with a conclusion that discusses the work that has been carried out and it provides directions for future work.

Main Principles

Each individual, or fly, corresponds to a 3D point that mimics a radioactive emitter, i.e. a stochastic simulation of annihilation events is performed to compute the fly's illumination pattern. For each annihilation event, a photon is emitted in a random direction. A second photon is then emitted in opposite direction. If both photons are detected by the scanner, the fly's illumination pattern is updated. The scanner properties (e.g. detector blocks and crystals positions) are modelled, and each fly is producing an adjustable number of annihilation events. Each fly keeps a record of its simulated LORs. Therefore the result of these simulations consist of a list, per fly, of pairs of detector identification numbers that correspond to LORs. These lists are aggregated to form the population total illumination pattern.

Initially, the flies' position is randomly generated in the volume within the scanner. Using genetic operations to optimise the position of radioactive emitters, the population of flies evolves so that the population total pattern matches measured data. The final population of flies corresponds to the tracer density in the patient, i.e. the reconstructed data. Note that cross-over operations are not used in this application. In our context, the result of such an operation may lead to meaningless results, e.g. in the case of cross-over between two flies of two distinct objects (the new fly will be wrongly located in between). Only mutation and immigration (i.e. a fly is created at a random position) are used.

Fitness Metrics

The fitness metrics corresponds to a distance measurement between the simulated data and the actual data given by the imaging system. City block distance provides a good compromise between accuracy and speed. Note that smaller the population's cost is, closer the simulated data is to the actual data.

In [START_REF] Bousquet | Fully three-dimensional tomographic evolutionary reconstruction in nuclear medicine[END_REF], we showed that, when we were addressing the SPECT problem, if we defined the fitness of a fly as the consistency of the image pattern it generates, with the actual images, it gave an important bias to the algorithm with a tendency of the smaller objects to disappear. This is why we then introduced marginal evaluation (F m (i)) to assess a given fly (i). It is based on the leaveone-out cross-validation method. We use a similar approach in PET:

F m (i) = dist (pop, input) -dist (pop -{i} , input) (1)
with F m (i) the marginal fitness of Fly i, dist (A, B) the city block distance between two tables A and B, pop is the set of LORs simulated by the whole population, input is the set of LORs extracted from the input data, and pop -{i} is the set of LORs simulated by the whole population without Fly i. The fitness of a given fly will only be positive when the global cost is lower (better) in presence rather than in the absence of this fly. We therefore used a fixed threshold to operate selection.

Thresholded Selection

At each iteration of the evolution loop, a fly has to be killed, and a fly may be used during the mutation. We saw in the previous section that the fly's fitness is its own contribution (positive or negative) with respect to the whole population.

We take advantage of this principle as follows: i) any "bad" fly (its fitness is negative) is a candidate for death, and ii) any "good" fly (its fitness is positive) is a candidate for mutation. When a fly is killed, its LORs are removed from the total set of simulated LORs. When a new fly is created, its LORs are added. This process needs to be fast to be able to decrease the number of bad flies and increase the number of good flies as much as possible.

Mitosis

To obtain accurate, high resolution images it is necessary to use large populations of flies. However, as processing time is roughly proportional to the number of flies to be processed, we choose a simple scheme that begins with a small population, then multiply the population along the algorithm execution using a mitosis process: each fly is duplicated. Newly created flies will have their own illumination pattern.

It is triggered whenever the number of flies with a negative fitness gets too low. In practice, at each step in the steady state process, one fly is chosen randomly and its fitness tested: a genetic operator will only be applied if the fly is bad. We launch the population mitosis every time 50 consecutive flies are found with a positive fitness.

Results

We have developed numerical phantom models to assess the reconstruction algorithm. To date, no scattering and no tissue attenuation have been considered. Whilst this is not physically correct, it allows us to test and validate our approach in the simplest cases. First, we present quantitative results in 2D-mode, both in low and high resolutions. Then, qualitative results in 3D-mode are presented using a complex object. For each test case, the initial population is 5,000 flies. When the current number of flies of the population is above a given threshold, e.g. 1 6 or 2 6 , the evolution loop is stopped whenever the number of flies with a negative fitness becomes too low, i.e. the stopping criteria is similar to the mitosis criteria. 70 LORs per fly have been simulated in the 2D test cases. The probability of LORs to be detected is much higher in the 3D case. Thus only 10 LORs per fly have been simulated in the 3D case.

2D-Mode

Test 1: large objects with different sizes and similar radioactivity concentration. The purpose of this test was to assess the ability of the algorithm to retrieve relatively large objects, whose sizes are different, but with the same radioactivity concentration. Fig. 1 shows the simulated set up. The phantom is made of two boxes (7 × 7 × 0.4 cm 3 and 10 × 10 × 0.4 cm 3 ) with the same radioactivity concentrations (∼ 930.000 counts/ml). The simulated PET system is made of a single ring of 72 linear blocks that include 8 crystals. To evaluate the results, a 512 × 512 pixel slice is produced (see Fig. 2). Note that the typical image size in PET is 128 × 128 pixels. The slices are post-filtered using a gaussian convolution kernel, then linearly rescaled between zero and one. Profiles in this reconstructed image are compared to corresponding profiles in the phantom data (see Fig. 3). For both boxes, the reconstructed data seems to be close to the input data. Full width at half maximum (FWHM) is also measured to quantify errors (see Table 1). These results show that our evolutionary scheme is able to accurately recover the width of our test objects.

Test 2: small objects with different sizes and radioactivity concentrations. This test case has been designed to assess the ability of our algorithm to detect small objects, and their relative radioactivity concentrations. etc.). A low resolution PET system has first been considered. It is made of a single ring of 72 linear blocks that include only 1 crystal.

To evaluate the results, a 512 × 512 pixel slice is produced once again (see Fig. 4(b)). The reconstructed data appears to be visually close to the input data. In particular, the size and concentration of cylinders are visually well preserved.

To estimate the diameter of each cylinder, horizontal profiles have been extracted so that they crossed the cylinders in their respective centre (see Fig. 5). As the lower profiles are symmetrically similar to the upper profiles, they are not plotted here. FWHM is measured once again to quantify errors (see Table 2). Let Object i be the cylinder whose concentration is C i in the phantom model and C ′ i in the reconstructed slice. Whilst the profiles in the reconstructed slice Table 2. FWHM estimated from Fig. 5 (using a low resolution PET system), and Fig. 6 (using a high resolution PET system).

FWHM from FWHM Relative difference FWHM Relative difference phantom model in Fig. 5 in Fig. 5 in Fig. 6 in Fig. (a) Phantom model. seem to match respective profiles in the phantom model, error measurements in FWHM are relatively high for the smallest cylinders (up to 35%). To investigate the influence of the reconstructed slice resolution with respect to the low spatial resolution of the PET system, the test case presented in the next section makes used of similar objects and a PET scanner with higher spatial resolution.

C ′ 1 C ′ 5 C ′ 3 C ′ 2 C ′
To assess the validity of the radioactivity concentration within cylinders, the average value at the centre of each cylinder has been measured in Fig. 4(b) (see Table 3). We compare the respective ratio of the different concentrations to the lower value (C ′ 1 ), so that we can compare the reconstructed values with the theoretical values. In theory, we should get 3 shows that the relative concentrations have been preserved in the reconstructed slice. However, the maximum relative error is about 16.50%.

C ′ 2 = 2C ′ 1 , C ′ 3 = 3C ′ 1 , etc. Table
Test 3: higher scanner resolution. The previous test case shows that our algorithm is able to retrieve the respective size of objects and their respective Table 3. Relative radioactivity concentration estimated from Fig. 4(b) (using a low resolution PET system), and Fig. 4(c) (using a high resolution PET system).

Relative

Relative error Relative Relative error concentration in Fig. 4(b) concentration in Fig. 4(c) Object in Fig. 4(b) (in %) in Fig. 4(c) concentration. However, relative errors can be as high as 35% for the FWHM and 16,5% for the concentration. In this test case, similar objects have been simulated. The size of crystals has been reduced so that their width matches the width of real crystals. A similar methodology is used to assess the results: i) Fig. 4(c) shows a 512 × 512 pixel slice that has been reconstructed, ii) profiles have been extracted (see Fig. 6), iii) FWHM estimated (see Table 2), and iv) the concentrations assessed (see Table 3). These results show that using a high resolution PET scanner reduced the maximum error in object size by a factor of 4. This is due to the improvement of the spatial resolution. On the other hand, errors in radioactivity concentration have not been reduced.

(in %) 1 C ′ 1 N/A C ′′ 1 N/A 2 2.13 × C ′ 1 6.5 2.17 × C ′′

3D-Mode

The last case has been performed in fully-3D, i.e. the PET imaging system is made of a stack of detector rings. For a coincidence event, the two photons of a LOR can be detected onto different rings. Only visual results are presented here. A complex 3D shape is used in this test. The simulation is performed using a polygon mesh (here we use the dragon model from The Stanford 3D Scanning Repository, http://graphics.stanford.edu/data/3Dscanrep/, last access 17 Jan 2010) that is uniformly filled with radio-tracers (see top row in Fig. 7). Then, LORs are recorded in fully-3D mode. Finally, we run our evolutionary reconstruction scheme. Note that the reconstruction algorithm is similar in both 2D and 3D modes. The only difference is the geometrical property of the simulated PET scanner. The bottom row in Fig. 7 presents the reconstructed dataset after volume rendering. One can visually distinguish the shape of the dragon from the population of flies.

Conclusion

It may occur that complex applications fuel fundamental technical developments. In the research presented here, we addressed a complex problem that had never been approached in the past using evolutionary computing, by transposing the Fly Algorithm technique originally developed in a stereovision context. We then faced several difficult issues which encouraged the development of new tools that can probably be used into other application fields in evolutionary computing. Using the 'marginal fitness' concept opened the way to using a simplified thresholded selection, which in turn allowed to introduce the mitosis operator that duplicates the population whenever the proportion of individuals with a negative contribution to the global fitness becomes too low, thus periodically reviving the efficiency of the classical operators (mutation and immigration).

Preliminary results on tests objects show the validity of this approach in both 2D and fully-3D modes. In particular, the size of objects, and their relative concentrations can be retrieved in the 2D mode. In fully-3D, complex shapes can be reconstructed.

To date, only true coincidence events have been considered. Further work will therefore include the use of more realistic input data (including random events and scattering), which will finally lead to implement the correction of scattering within our algorithm. A comparison study against ML-EM and/or OS-EM methods will also need to be conducted.

Fig. 1 .

 1 Fig. 1. Simulated PET System: a single ring of 72 linear blocks that include 8 crystals; two boxes (7 × 7 × 0.4 cm 3 and 10 × 10 × 0.4 cm 3 ) with the same radioactivity concentrations (∼ 930.000 counts/ml).

Fig. 2 .

 2 Fig. 2. Slices (512 × 512 pixels) through the cubes.

Fig. 4 Fig. 3 .

 43 Fig. 3. Profiles extracted from Fig. 2.

4 (

 4

Fig. 4 .

 4 Fig. 4. Slices (512 × 512 pixels) through the cylinders.

Fig. 5 .

 5 Fig. 5. Profiles extracted from Fig. 4(b).

  2 nd line.

Fig. 6 .

 6 Fig. 6. Profiles extracted from Fig. 4(c).

Fig. 7 .

 7 Fig. 7. Top row: simulated object; bottom row: volume rendering of the reconstruction.

Table 1 .

 1 FWHM estimated from Fig.3

	FWHM from phantom model FWHM from slice Relative difference
	(in mm)	(in mm)	(in %)
	71	72	1
	99	99	0

⋆ member of Fondation Digiteo (http://www.digiteo.fr/). Now with University of California San Diego, Rebecca and John Moores Cancer Center, Radiation Oncology, CA, USA.


Franck P Vidal
email: franck.p.vidal@gmail.com

Evelyne Lutton
email: evelyne.lutton@inria.fr

Jean Louchet
email: jean.louchet@gmail.com

Jean-Marie Rocchisani
email: jean-marie.rocchisani@avc.aphp.fr

Threshold selection, mitosis and dual mutation in cooperative co-evolution: application to medical 3D tomography

We present and analyse the behaviour of specialised operators designed for cooperative coevolution strategy in the framework of 3D tomographic PET reconstruction. The basis is a simple cooperative co-evolution scheme (the "fly algorithm"), which embeds the searched solution in the whole population, letting each individual be only a part of the solution. An individual, or fly, is a 3D point that emits positrons. Using a cooperative co-evolution scheme to optimize the position of positrons, the population of flies evolves so that the data estimated from flies matches measured data. The final population approximates the radioactivity concentration. In this paper, three operators are proposed, threshold selection, mitosis and dual mutation, and their impact on the algorithm efficiency is experimentally analysed on a controlled test-case. Their extension to other cooperative co-evolution schemes is discussed.

Introduction

Evolutionary algorithms have been proven efficient to solve the inverse problem of 3D data reconstruction in tomography [START_REF] Bousquet | Fully three-dimensional tomographic evolutionary reconstruction in nuclear medicine[END_REF], and particularly of positron emission tomography (PET) reconstruction in nuclear medicine [START_REF] Vidal | Artificial evolution for 3D PET reconstruction[END_REF][START_REF] Vidal | PET reconstruction using a cooperative coevolution strategy in LOR space[END_REF][START_REF] Vidal | New genetic operators in the fly algorithm: application to medical PET image reconstruction[END_REF].

In PET, a positron emitter is used as radionuclide for labelling. Positrons generally lead to an annihilation reaction, that emits two photons of 511 keV in opposite directions. This radiation is detected in coincidence, i.e. using the difference in arrival times of the detected photons of each pair, and considering that each annihilation produces two photons emitted in exactly opposite directions. The line joining the detectors that have been activated for a given pair of photons is called "line of response" (LOR). An overview of reconstruction methods in nuclear medicine can be found in [START_REF] Zeng | Image reconstruction -a tutorial[END_REF].

In previous work, we showed that a cooperative coevolution strategy (or Parisian evolution) called "flies algorithm" [START_REF] Louchet | Stereo analysis using individual evolution strategy[END_REF] could be used in Single-Photon Emission Computed Tomography (SPECT) reconstruction [START_REF] Bousquet | Fully three-dimensional tomographic evolutionary reconstruction in nuclear medicine[END_REF], and also in PET reconstruction in 2D-mode [START_REF] Vidal | Artificial evolution for 3D PET reconstruction[END_REF][START_REF] Vidal | PET reconstruction using a cooperative coevolution strategy in LOR space[END_REF], and in Fully-3D-mode [START_REF] Vidal | New genetic operators in the fly algorithm: application to medical PET image reconstruction[END_REF]. The marginal fitness was used to propose new operators, i) the threshold selection, and ii) the mitosis, but no performance analysis of these operators has been performed so far.

This paper addresses this deficiency and it analyses the impact of each operator on the performance of a PET reconstruction algorithm on a controlled test-case. A new operator (namely the dual mutation) and a pre-initialisation of the flies' position using back-projection are also described and analysed. Standard PET reconstruction algorithms are reviewed in Section 2. It is followed by an overview of the fly algorithm for PET reconstruction. The three operators that control a varying population size scheme are presented in Section 4, as well as an alternate initialisation process. Experimental setup and analysis are given in Section 5, before presenting some conclusions and future work in Section 6.

Standard PET reconstruction algorithms

Tomography reconstruction algorithms can be divided into two main categories.

On the one hand, there are analytical methods. These are based on a continuous modelling and the reconstruction consists in the inversion of measurement equations, such as the well known Filtered Back-Projection (FBP). This method is now rarely used due to strong artefacts in the reconstructed data (see Fig. 5) and also because the correction of imaging physics effects need to be undertaken before the reconstruction, leading to a systematic positive bias in the reconstructed volume.

On the other hand, there are iterative methods. This class of methods can be split into two kinds. Algebraic methods are used in X-ray Computed Tomography (CT); statistical methods are used in nuclear medicine for both SPECT and PET [START_REF] Vandenberghe | Iterative reconstruction algorithms in nuclear medicine[END_REF]. They take into account noise, and the correction of imaging physics can be applied during the reconstruction in the iterative steps. Iterative methods are relatively easy to model. In practice, the volume is usually discretised into voxels. Each voxel intensity is treated as an unknown. A system of linear equations is defined according to the imaging geometry and physics: p = R f , with f the volume to recover, p the measured data, R the system model. Imaging physics, such as non-uniform attenuation, scatter, etc. can be modelled in R, whereas they are difficult to handle in an analytic algorithm. The system of equations is finally solved using the iterative algorithm.

There are different ways to implement these iterative methods. The main differences are about the computation of the projections, the physics corrections (scattering, random, attenuation, etc.) are applied, and how the error corrections are applied in the estimated projections.

The Maximum Likelihood -Expectation Maximisation (ML-EM) (or 'EM") is a common algorithm in SPECT and PET. It assumes Poisson noise is present in the projection data. ML-EM does not produce artefacts seen in FBP reconstructions, and it has a better signal-to-noise ratio in region of low concentration. However, the algorithm converges slowly.

The Ordered Subset -Expectation Maximization (OS-EM) has been proposed to speed-up convergence of the EM algorithm. Its principle is to reduce the number of projections used at each iteration of the EM algorithm. Projections are grouped in K sub-groups. The projections of a sub-group are uniformly distributed around the volume to reconstruct.

PET reconstruction using the fly algorithm

The fly algorithm for tomography reconstruction follows the iterative paradigm. The steps of the iterative method can be described as follows:

1. Each individual, or fly, corresponds to a 3D point. Initially, the flies' position is randomly generated in the volume within the scanner. The population of flies corresponds to the tracer density in the patient. 2. To produce estimated projection data, each fly mimics a radioactive emitter, i.e. a stochastic simulation of annihilation events is performed. For each annihilation event, a photon is emitted in a random direction. A second photon is then emitted in the opposite direction. If both photons are detected by the scanner, the corresponding LOR is recorded. The scanner properties (e.g. detector blocks and crystals positions) are modelled, and each fly is producing an adjustable number of annihilation events. 3. The optimisation is performed using genetic operations. The fitness function used during the selection operation takes into account the comparison between the estimated projections and the measured projections. 4. Using genetic operations to optimise the position of radioactive emitters, the population of flies evolves so that the population total pattern matches measured data. 5. Instead of a "generational" evolutionary strategy, in which at each loop every individual (fly) will be eliminated and replaced with a new fly, we chose a "steady state" evolutionary strategy.

Note that in classical evolutionary approaches, each individual in the population is a potential solution; in the Fly approach, a subset of the evolving population itself is the representation of the solution. After convergence, the "good" flies (see Section 4.1) are then extracted to form the reconstructed volume.

Varying population size scheme in a cooperative co-evolution algorithm

Cooperative co-evolution strategies rely on a "social" formulation of the optimisation problem, where individuals collaborate or compete in order to collectively build a solution. The fly algorithm is a mono-population strategy (Parisian approach): all flies contribute independently and collectively to build the solution.

In [START_REF] Barriere | Experimental analysis of a variable size mono-population cooperative-coevolution strategy[END_REF] a variable sized population Parisian GP strategy has been successfully used on a cooperative co-evolution, based on adaptive population deflating and inflating schemes. We test in this paper an "inflating-only" strategy, the mitosis, described below, to gradually increase the precision of the reconstructed data.

Marginal fitness

In this application, the similarities or discrepancies between the estimated projection data and the measured projection data provided by the imaging system have to be assessed. We chose City Block distance (also known as Manhattan distance) as the fitness metrics to measure the distance between two LOR sets. It provides a good compromise between speed and accuracy. Eq. 1 provides the global fitness, i.e. the population's cost:

dist(LOR m , LOR e ) = M i M j |LOR m (i, j) -LOR e (i, j)| (1)
with dist(LOR m , LOR e ) the City Block distance between LOR m and LOR e , the set of LORs for the measured data and the estimated data respectively, LOR(i, j) is the number of counts of a LOR between the photon detectors i and j, M is the total number of photon detectors within the imaging system. LOR sets are efficiently implemented using triangular sparse matrices to reduce the amount of memory needed to store the data. The smaller global fitness is, the closer the simulated data will be to the measured data.

In [START_REF] Bousquet | Fully three-dimensional tomographic evolutionary reconstruction in nuclear medicine[END_REF], we showed that, when we were addressing the SPECT problem, if we defined the fitness of a fly as the consistency of the image pattern it generates, with the actual images, it gave an important bias to the algorithm with a tendency of the smaller objects to disappear. To address this, we introduced marginal evaluation to assess a given fly. We use a similar approach in PET:

F m (x) = dist (LOR e -{LOR x } , LOR m ) -dist (LOR e , LOR m ) (2)
with F m (x) the marginal fitness of Fly x, and LOR e -{LOR x } is the set of LORs simulated by the whole population without Fly x. In practice, each fly needs to keep a record of its simulated LORs. The fitness of a given fly will only be positive when the global cost is lower (better) in presence rather than in the absence of this fly.

Threshold selection

The fly to be killed is randomly chosen by the "selection" operator, with a bias towards killing "bad" individuals. On the other hand, if the new fly is to be created by mutation of another fly, this fly is randomly chosen by the "selection" operator, with a bias towards reproducing "good" individuals. Classical selection operators are ranking, roulette wheel and tournament [START_REF]Evolutionary Computation 1: Basic Algorithms and Operators[END_REF]. In our algorithm, as each fly's fitness is the value of its (negative or positive) contribution to the quality of the whole population, we managed to simplify and speed up the selection process by using a fixed fitness threshold. Any "bad" fly (its fitness is negative) is a candidate for death, and any "good" fly (its fitness is positive) is a candidate for mutation.

Mitosis

When the number of flies with a negative fitness decreases, the threshold selection fails to provide flies to be killed in an acceptable time. It also means that the reconstruction is optimum at the current resolution. If the resolution is acceptable, i.e. there are enough flies to approximate the radio-tracer concentration, then the algorithm can stop and the reconstructed volume is extracted using flies with a positive fitness. If not, a mitosis operator is triggered to gradually increase the population size. Each fly is split into two new flies to double the population size. One of the two flies is then mutated.

Dual mutation

To optimise the flies' position, our algorithm takes advantage of a mutation operator. When a new fly (b) is created by mutation of an old "good" fly (a), the position of Fly b is first initialised to the same position as Fly a. The new fly is then stochastically translated in any direction, and LORs are randomly generated from that fly. The length of the translation vector is a random variable that follows a Gaussian law whose mutation variance is σ 2 . It needs first to be set to a large value to better explore the search space. However, a constant large mutation variance will lead to blurred reconstructed volumes. σ has therefore to be gradually reduced.

The use of adaptive mutations in evolutionary algorithms is an ancient idea, directly inspired by natural adaptive phenomena, e.g. mutations simulated by stress [START_REF] Rosenberg | Evolving responsively: adaptive mutation[END_REF]. In artificial evolution, various adaptive schemes have been considered for mutation [START_REF] Eiben | Parameter control in evolutionary algorithms[END_REF], depending of the parameter to be adapted (standard deviation, σ [START_REF] Bäck | An overview of evolutionary algorithms for parameter optimization[END_REF], or mutation law [START_REF] Chellapilla | Combining mutation operators in evolutionary programming[END_REF] for continuous mutation, mutation probability for discrete mutations [START_REF] Ochoa | Setting the mutation rate: Scope and limitations of the 1/L heuristic[END_REF]). Regarding the adaptation of σ, one can distinguish several strategies :

σ is directly adapted to local measurements, like fitness [START_REF] Lis | Application of evolutionary algorithms with adaptive mutation to the identification of induction motor parameters at standstill[END_REF] or local regularity [START_REF] Lutton | Pointwise regularity of fitness landscapes and the performance of a simple ES[END_REF], -σ is tuned depending on some success measurement: in this category fall the famous 1/5 th rule proposed by Schewefel [START_REF] Schwefel | Numerical optimization of computer models[END_REF][START_REF] Beyer | Evolution strategies -a comprehensive introduction[END_REF], -σ is subject to an adaptive pressure itself, it is self-adapted [START_REF] Bäck | Self-adaptation in genetic algorithms[END_REF]: σ is considered as an additional parameter in the genome, and a log-normal Gaussian law is used to control the "mutation over the mutation".

These techniques have been proven efficient in various cases, depending on the fitness function and the genetic engine. It has however to be noticed that the sophistication of a mutation operator has a computational cost, and that some very rough schemes may perform better due to their capability to rapidly test numerous sample points [START_REF] Collet | Issues on the optimisation of evolutionary algorithm code[END_REF].

Concurrent testing with various subpopulation has been also considered for mutation law adaptation [START_REF] Chellapilla | Combining mutation operators in evolutionary programming[END_REF].

Here we propose an adaptive mutation scheme based on the concurrent testing of two alternative σ values (σ low and σ high , with k σ low = σ high ). The update rule is multiplicative as for the 1/5 th rule. If σ high gives the best results during the previous period, then both mutation variances are increased by a predefined factor (pf , with pf > 1). If σ low gives better results, then the variances are multiplied by 1 pf . The major advantage of this dual mutation scheme is to provide a fully automatic method to adapt the mutation variance, whilst keeping the administration cost of the algorithm relatively light. Additionnally, this scheme does not need to make any assumption on the ideal success rate of the mutation as in the 1/5 th rule. In practice, the global fitness is recorded after each mutation. The cumulative difference of the global fitness (∆(σ)) before and after mutations is computed to determine which σ value provides the best performance over a given period of time. To prevent oscillation of σ values, a criteria can be added to avoid changes when both σ high and σ low provide relatively similar results, e.g. when the absolute difference between ∆(σ low ) and ∆(σ high ), relative to the current global fitness, is below a given threshold (t mut ).

Initialisation of flies on LORs

Iterative reconstruction methods generally make use of a constant volume as an initial estimate of the volume (see Fig. 5(a)).

However, to speed-up the reconstruction process, a volume is first reconstructed using a fast analytical algorithm, the simple back-projection, that we implemented on the graphics card using OpenGL. The algorithm consists in back-projecting each LOR into the volume space. Pixels along the path of a LOR are updated uniformly. This operation is fast and provides the evolutionary algorithm with an initial guess of the volume (see Fig. 5(c)). For each voxel of the initial estimate, a given number of flies is assigned depending on the voxel intensity (see Fig. 5(b)).

Results

The validity of the reconstruction method has been addressed in [START_REF] Vidal | New genetic operators in the fly algorithm: application to medical PET image reconstruction[END_REF]. In this paper, we focus on the evaluation of the performance of the new genetic operators. For each test case, 750000 new individuals have been created. For each tested configuration, the reconstruction has been repeated 20 times, and the final global fitness was recorded. For every test, unless specified, the dual variance and the threshold selection operators have been enabled. Results are presented using box plots (also called box-and-whisker diagrams).

Experimental setup

Here, a single ring PET system is considered. Its radius is about 430mm. The ring is made of 72 linear blocks that include 8 crystals each. The width of a crystal is about 4.5mm. Fig. 1(a) shows the reference image. It includes nine cylinders having two different radii (1 cm and 2.5 cm) and five different radioactivity concentrations (C 1 = 114, 590 count/ml,

C 2 = 2C 1 , C 3 = 3C 1 , etc.) C1 C5 C3 C2 C4
(a) Reference image.

C ′′ 1 C ′′ 5 C ′′ 3 C ′′ 2 C ′′ 4 (b) Reconstructed image.

Threshold selection

The size of the population is fixed (160000 flies), i.e. no mitosis has been triggered. The performance of the threshold selection and the tournament selection are presented in Fig. 2. Both operators provide similar performance. The threshold selection is then preferred because of the additional information that it brings: enable mitosis, and provide a convergence criteria at a given resolution (i.e. for a given size of population).

Mitosis

Two variables have to be assessed at the end of the reconstruction: the current size of the population, and the global fitness. The larger the final population, the better the image resolution can be obtained. The smaller the global fitness, the closer the estimated data to the measured data is. Fig. 3(a) shows the average number of flies in the final population depending on the initial size of the population (625, 2500, 10000, 40000, 80000, and 160000 flies). When the size of the population is 160000 flies, no mitosis has been triggered. with the smallest initial population size. Then, the reconstruction converges much faster when the initial population is small (smaller global fitness and bigger final population size). These results validate the efficiency of the mitosis operator.

Dual mutation

The initial σ low value in this test is 35mm, pf is equal to 3 √ 2, and σ high is equal to prevent any change of σ values, leading to unsatisfactory results. Fig. 4(b) shows the global fitness obtained i) using a constant variance (see ( 1) and ( 2)), or ii) enabling dual mutation operator (see ( 3) and ( 4)). The best results are observed using the dual mutation operator with a very low t mut value. the imaging system. In the latter case, the position of flies is initialised using the simple back projection. Fig. 3 shows the performance of both strategies when the mitosis operator is enabled.

Initialisation of flies on LORs

When the initial size of the population is relatively large, the algorithm converges much faster using this initialisation step. This is not the case when the initial size of the population is relatively small. It may be due to the fact that the algorithm converges fast enough when only a few flies are used. When the initial number of flies is slightly higher, the reconstruction converges faster when the position of flies are initialised using the back projection.

Conclusion and futher works

We have presented new operators in cooperative co-evolution and validated their efficiency using a controlled test-case in PET reconstruction. Both the threshold selection, mitosis and dual mutation operators have shown their usefulness and ability. Experimental statistics show that the threshold selection perform as well as the tournament selection, but it has the great advantage of bringing a convergence criterion related to the current resolution. Additionally, it allows to trigger an automatic mitosis, i.e. doubling the population size, to improve the resolution. Best performance, both in term of final resolution and convergence, are obtained using small initial population size. The dual mutation operator provides an adaptive mutation variance that has proven to be better than using fixed mutation variances.

Such operators can be used in any other cooperative co-evolution schemes as soon as a marginal fitness can be considered as beneficial, that obviously depends on the computation cost of the marginal fitness. For instance, threshold selection, mitosis and dual selection will be considered as further work for the original fly algorithm on a stereo-vision application ( [START_REF] Louchet | Stereo analysis using individual evolution strategy[END_REF]). The marginal fitness will also be considered for developing a "deflating operator". This additional mechanism for controlling the population size may be interesting in the case of applications whose resolution does not depend on the size of the final population. Further work will also include the correction of photon attenuation and Compton scattering, and a concurrent study with the OS-EM algorithm.

Fig. 1 .

 1 Fig. 1. Slices (512 × 512 pixels) through the cylinders.

Fig. 2 .

 2 Fig. 2. Performance of the threshold selection and of the tournament selection.

Fig. 3 (Fig. 3 .

 33 Fig.3. Performance of the mitosis operator using variable initial population sizes.

3 √Fig. 4 .

 34 Fig. 4. Performance of the dual variance operator.

Fig. 5 (

 5 Fig. 5(a) and Fig. 5(b) show two possible initial estimates of the radio-active concentration. In the first case, flies are uniformly located within the space in

  (a) Uniform distribution of the flies. (b) Initialisation of the flies' position using Fig. 5(c) (c) Image reconstructed using the simple back projection algorithm.

Fig. 5 .

 5 Fig. 5. Initial estimates of the reconstructed image.


Salma Mesmoudi
email: salma.mesmoudi@grignon.inra.fr

Nathalie Perrot
email: nathalie.perrot@grignon.inra.fr

Romain Reuillon
email: romain.reuillon@iscpif.fr

Paul Bourgine
email: paul.bourgine@polytechnique.fr

Evelyne Lutton
email: evelyne.lutton@inria.fr

OPTIMAL VIABLE PATH SEARCH FOR A CHEESE RIPENING PROCESS USING A MULTI-OBJECTIVE EA

Keywords: Multiobjective Evolutionary aglorithm, viability modeling, optimal path search, indirect encoding, agri-food process modeling, cheese ripening

Viability theory is a very attractive theoretical approach for the modeling of complex dynamical systems. However, its scope of application is limited due to the high computational power it necessitates. Evolutionary computation is a convenient way to address some issues related to this theory. In this paper, we present a multiobjective evolutionary approach to address the optimisation problem related to the computation of optimal command profiles of a complex process. The application we address here is a real size problem from dairy industry, the modeling of a Camembert cheese ripening process. We have developed a parallel implementation of a multiobjective EA that has produced a Pareto front of optimal control profiles (or trajectories), with respect to four objectives. The Pareto front was then analysed by an expert who selected a interesting compromise, yielding a new control profile that seems promising for industrial applications.

Introduction

In this paper, we consider a multi-objective optimisation problem related to the modeling of an agrifood industrial process, that is the camembert cheese ripening process. Experimental and analytical approaches allow to build a set of models that reflect the numerous interactions that take place at different scales, from microscopic to macroscopic levels. These models are combined into a single complex one, that is exploited via as a so-called viability analysis.

Viability theory considers the problem of maintaining a system under a set of constraints [START_REF] Aubin | Viability Theory[END_REF]. To solve viability problems one generally builds a viability kernel, which is the set of initial states from which there exists a trajectory that remains in the set of constraints. This theory is very convenient to model complex systems, but current numerical applications are limited due to the high computational power it necessitates. A question currently addressed by the community is the computation of a viability kernel, and various approaches now exists [START_REF] Bokanowski | An anti-diffusive scheme for viability problems[END_REF].

We deal here with another critical problem related to the application of this theory, which is how to practically use a viability kernel once it has been computed. A question that needs to be addressed is the one of finding an optimal path within the viability kernel. By optimal path we mean a set of control parameters and commands that guides the dynamical system toward a desired target state (that is, in the case of a cheese, a desired target weight and quality). We will see in the sequel that it can be formulated as a multiobjective optimisation problem.

The paper is organised as follows. Section 2 describes some basic principles of a viability analysis applied to cheese ripening modeling, and describes the optimisation problem that has to be solved. Then the OpenMOLE platform on which algorithms have been developed, is described in section 3. Section 4 presents the details of the algorithm that has been developed: it is based on an indirect encoding of a path in the viability graph. Experimental setup and results are given in section 5. An analysis of the optimal trajectory is then made by a cheese-ripening expert in section 5.3, before concluding (section 6). cheese ripening process Industrial cheesemaking of camembert is based on the use of pasteurized milk, and on a ripening process conducted in controlled chambers, in order to obtain a product at a given (savety and quality) target.

A camembert type cheese is made by inoculating milk with lactic acid bacteria (Flora Danica lyophilisate, CHN11, Chr Hansen, Arpajon, France), and other specific microorganisms like Kluyveromyces marxianus (GMPA collection, 448), Geotrichum candidum (Degussa,D), Penicillium camemberti (Degussa,R), and Brevibacterium aurantiacum (ATCC9175) as ripening flora. After coagulation and cutting, it is then shaped in moulds. After drained (around 24h), fresh cheese is obtained and transfered in the ripening chamber. At this stage, each cheese weights around 300g. The cheeses are then left for ripening for at least 3 weeks, in order to obtain to obtain (1) a limited mass loss and at the same time (2) a given quality of coat of the rind and a creamy texture, due to an optimum microorganisms behavior.

For modeling purpose, we consider here a critical phenomenon for industrials which is the cheese mass loss during the ripening process. It is linked in a complex way to evaporation and to carbon consumption due respiration of micro-organisms in the ripening chamber [START_REF] Helias | Modeling of camembert-type cheese mass loss in a ripening chamber. Main biological and physical phenomena[END_REF], for instance:

• Low relative humidity and high temperature on the cheese surface increases evaporation.

• Cheese surface temperature decreases when evaporation occurs.

• Respiration increases the cheese surface temperature as heat is produced during the substrate degradation.

The state variables that have then been chosen to build the viability kernel are:

• the cheese mass, with a range from 250g to 310g with steps of 1g,

• the cheese temperature from 8

• C to 16 • C (step 1 • C)
• the respiration rate of micro-organisms from 0 to 50 g/m 2 /day (step 1g/m 2 /day).

Additionnally, the control variables are

• the ripening room temperature (from 8

• C to 16 • C, step 1 • C),
• the relative humidity (from 84% to 98%, step 2%).

The viability kernel is defined by the following constraints:

• The ripened cheese should have at the end of the process a mass between 250-270g, a temperature between 8 • -10 • C and a respiration between 23 and 50g/m 2 /day.

• The (CO 2 , O 2 ) gas rate evolution should follow a specified profile along the process, in order to fully exploit the micro-organisms capabilities. This constraint is set by experts. The respiration rate should begin at level 0 on day 1 (microbial growth latency), reaches a maximum between day 3 and day 8 and decreases slowly during the last days of ripening (see [START_REF] Sicard | The viability theory to control complex food processes[END_REF] for more details).

On the basis of this model, a viability kernel has been calculated [START_REF] Sicard | The viability theory to control complex food processes[END_REF]. Coupled to this calculus, a geometric analysis of the shape of the viability kernel, and of all the viable paths to the target (also named the viability tube), provides useful informations about the robustness and uncertainties of the system. This analysis is based on the computation of the boundaries of the viability tube, and for each point of the tube, on its distance to the nearest boudary. Optimal algorithms for the Euclidean distance transform (EDT) in arbitrary dimension have been developed for morphological mathematics and image analysis purpose. [START_REF] Mesmoudi | Geometric analysis of a capture basin. application to cheese rippening process[END_REF] adapted one of these algorithms [START_REF] Coeurjolly | Optimal separable algorithms to compute the reverse euclidean distance transformation and discrete medial axis in arbitrary dimension[END_REF] for the analysis of viability tube data.

The aim is to find an optimal strategy, i.e. an optimal and real path, among each possible viable path to the target. However each trajectory in the viability tube can be evaluated according to four goals : loss mass minimisation, number of control changes, trajectory robustness and optimal breath evaluated by the time to reach the maximal respiration rate. Convenient algorithms for that purpose are thus multiobjective optimisation ones.

Multi-objective optimisation problems are often NP-hard, complex and CPU time consuming. Exact methods can be used to find the exact Pareto front (or a subset of the front), but usually it is impossible to compute exact solutions for large problems, as they are time and memory consuming. For instance, for the cheese ripening case, there exists 10 24 possible trajectories in the viability kernel we defined.

We rely in this work on the use of multi-objective Evolutionary Algorithm (EA) to approximate Pareto fronts within a reasonable computation time. Additionnally, as the cheese mass loss modeling is actually a large sized problem, the algorithms have been implemented on a parallel and distributed computation grid.

Our problem corresponds to a search for an optimal path in a graph, where each node represents the state of the cheese (mass and breath) and each edge represents a control (relative humidity and chamber temperature). The graph is organised as a succession of slices, each representing a day, it is therefore not complete.

Parallel Implementation

In this paper, we deal with a grid of computers (multicomputers architecture), a multiple-population (or multiple-deme) structure is more convenient. These algorithms maintain several subpopulations that occasionally exchange individuals (during a migration step). Additionnally, we choose to use the multiobjective approach proposed by [START_REF] Horn | A niched pareto genetic algorithm for multiobjective optimization[END_REF] and [START_REF] Deb | A fast elitist non dominated sorting genetic algorithm for multi-objective optimization: Nsga-ii[END_REF]: Nondominated Sorting Genetic Algorithm (NSGA) and NSGAII.

Our implementation is based on the OpenMOLE platform of "workflow" [START_REF] Reuillon | Declarative task delegation in openmole[END_REF]. The development of OpenMOLE was organized under the form of a community of free software programs that is today very active. OpenMOLE is a framework providing distributed computing facilities. It takes advantage of generic interface of JSAGA (http://grid.in2p3.fr/jsaga/) and provides on top of that. OpenMOLE has been designed to work out of the box on the user desktop with the idea of completely hiding the fact that computation may be carried out on distributed environments.

Algorithm

The NSGA-II algorithm is based on a ranking that shares the population into several classes at each gen-eration (see [START_REF] Deb | A fast elitist non dominated sorting genetic algorithm for multi-objective optimization: Nsga-ii[END_REF] for more details).

Search space, solutions encoding

The Camembert ripening process can be described via a graph, where each node is a viable state, i.e. a point in a space of dimension 3: mass(M) (in g), breathing (R) (in g/m 2 /day) and time (i) (in days of ripening). This graph is oriented with respect to time. It is incomplete since every state cannot be connected to all others states. Usually, for every viable state there exists k outgoing edges, or controls, such as 1 ≤ k ≤ 74. Every viable state E M i ,R i ,i such as 2 ≤ i ≤ N is at least connected by two oriented edges:

E M i-1 ,R i-1 ,i-1 → E M i ,R i ,i → E M i+1 ,R i+1 ,i+1
to each edge (→) corresponds a control pair C i k : that represents the relative humidity and the ripening room temperature to be applied to obtain the connected node (objective state), for instance

E M i+1 ,R i+1 ,i+1 from E M i ,R i ,i .
A first question to address now is how to encode a viable path in a genome and respect the various validity constraints. A direct encoding is for instance an ordered list of states, that represent the successive states of a camembert during its ripening. Each state change, that corresponds to an edge of the viability graph, is characterised by its controls C i k . This representation has however an important drawback, as it becomes difficult to maintain a set of valid solutions after the application of genetic operators.

We thus preferred an indirect encoding, i.e. to encode the way a graph is built instead of encoding the graph itself. Any viable path can be specifyed, by its initial state and by a set of successive controls. But, as the the number of outgoing controls/edges may be different for each state, it has been necessary to built an ad-hoc addressing that is valid whatever the state.

Let C i = C i 1 ,C i 2 , . . . ,C i k be the set of possible outgoing control pairs for the state E M i ,R i ,i (k is variable and depends on the current state). To specify a given edge it is thus enough to provide an index that identify the edge in the set C i .

A robust way to do this, while ensuring a valid edge is uniquely identified at each state, is to use an edge addressing based on real numbers, as follows. Let S = {p(C 1 ), p(C 2 ), .., p(C N )} be a set of reals of [0, 1] where each p(C i ) represents a proportion, or percentage that is interpreted with respect to the number of outgoing edges of the state E i . For instance, if in a state E i has 10 outgoing edges, then if p(C i ) belongs to [0, 1 10 ] it represents the first edge, the second one if it belongs to ] 1 10 , 2 10 ], so on until ] n-1 10 , n 10 ] for the last edge. The edges are ordered using the values of the associated control pair (humidity and temperature).

Fitness functions, pareto optimality

An optimal ripening path can be defined with respect to four objectives: 1. mass loss, to be minimised 2. number of control changes, to be minimised, 3. trajectory robustness, to be maximised, 4. number of days to get the maximum of respiration, to be minimized. This problem is a typical multiobjective optimisation, usually defined as follows.

"min/max"z = f (x) = ( f 1 (x), f 2 (x), ..., f m (x)) ∈ R m with x = (x 1 , x 2 , ..., x n ) ∈ X, an n -dimensional vec-
tor, and X the search space. f i are partial evaluations of the solution, usually corresponding to contradictory aims. The optimal solutions correspond thus to a set of compromise between the m various partial evaluations. In other words, the Pareto optimal set X * is made of all non-dominated points, i.e. points for which it is impossible to improve any objective without simultaneously worsening another:

X * = {x * ∈ X| ∃x ∈ X, f (x) ≤ f (x * )}, where f (x) ≤ f (y) ⇐⇒ ∀i ∈ 1..m, f i (x) ≤ f i (y)
and ∃ j ∈ 1..m, f j (x) < f j (y)

Genetic engine

• Selection and elitism: Among the various selection operators that have been proposed in the literature (see for instance [START_REF] Goldberg | Genetic algorithm in search, optimization and machine learning[END_REF])), we choose to use a tournament selection, due to its parsimonious mechanism that only considers a small random subsample of the population. Additionnally, comparisons are based on both domination and crowding criteria [START_REF] Tsai | A modied multiplesearching method to genetic algorithms for solving traveling salesman problem[END_REF]. • Crossover: In the case of a direct coding, i.e.

when the path is encoded as a set of successive states, it becomes necessary to repair invalid genomes after crossover. There exists various solutions in the literature to cope with orderbased encodings, like cycle crossover [START_REF] Oliver | A study of permutation crossover operators on the traveling salesman problem[END_REF], partially matched crossover (Goldberg, 1989) and order crossover [START_REF] Goldberg | Genetic algorithm in search, optimization and machine learning[END_REF], [START_REF] Davis | Applying adapting algorithms to epistatic domains[END_REF]. However, the indirect encoding we propose in section 4.1 allows to directly use simple crossovers such as the uniform crossover (UX) [START_REF] Syswerda | Uniform crossover in genetic algorithms[END_REF] and arithmetic crossover (AX) [START_REF] Eiben | Introduction to Evolutionary Algorithms[END_REF]. • Mutation: In this work we used a basic random mutation.

Results

Parameters setting

The parameter values were set heuristically. Additionally, as the computational load is high (around 8 hours using a computionnal grid) we did not perform a comprehensive study on the influence of different parameter settings, and it is possible that a careful fine-tuning of some values could bring slight improvements to the achieved results. For all experiments, the parameter setting is given in table 1.

Numerical results

The Pareto front that has been estimated is made of 142937 trajectories, however, with respect to the objectives, they correspond to only 329 combinations (see figure 2). A principal component analysis (PCA) was performed on the Pareto front obtained using our implementation. The PCA was based on the following variables: mass loss minimisation, number of control changes, trajectory robustness and optimal breath evaluated by the time to reach the maximal respiration rate. The first two eigenvectors represent 45.35% and 28% of the total variance, respectively. The variable projection and the distribution of the Pareto front solutions are represented in Figure 2. we can notice that the objective space is diversified: several solutions are associated to each objectif.

Analysis

This Pareto front has been analysed by a cheese ripening expert, who selected an efficient optimised trajec- The optimised trajectory (PT), has been compared to a viable trajectory (TVA) computed in a 8-day viability kernel and already applied a pilot [START_REF] Sicard | Development of a viability approach for reverse engineering in complex food processes: Application to a camembert cheese ripening process[END_REF], and to a standard one (SRT) running in 12day and used in dairy industry. This TP trajectory differs from the TVA trajectory and from the classical one. The relative humidity is not constant like in TVA (94%) or in the standard SRT trajectory (92%) (table 2). However, like in TVA (table 2) the temperature control varies whereas in the standard trajectory it remains constant. To analyse the consequences of these PT trajectory control changes, we compare its cheese mass loss evolution and respiration rate to those of TVA and standard trajectories.

Figure 3 shows the expected mass loss during the PT trajectory compared to the mass loss during the TVA and SRT trajectories. The quantities that are displayed are computed for PT, and measured during an experimental real simulation in a ripening chamber for TVA and SRT. The mass loss is 0.013 kg for the PT trajectory, while it is of 0.034 kg for the TVA tra- jectory and of 0.054 kg for the standard ripening. This result is significative, as a minimisation of mass loss is a very challenging issue for the dairy industry. This small mass loss can be explained by the high values of relative humidity of the PT trajectory. The cheese mass at the end of the robust ripening remains within the desired target range (0.25 kg to 0.27 kg).

Nevertheless, the mass loss is not enough to assess the quality of a cheese ripening process. Mass loss is a very sensitive parameter (a mass can be lost in a dried atmosphere in a few days only), but it is not the only one. If the other phenomena are not correctly controlled, ripening may not be satisfying. The microbial activities of optimized (PT), viable (TVA) and standard (SRT) ripening processes were also compared, and the results are presented in Figure 4. The respiration rate is the only kinetic that was checked. As shown, the respiration rate starts at 0, reaches a maximum of over 40 g/m 2 /day, and then slowly decreases until the day the cheese is wrapped. The maximum respiration rate begins two days earlier in the PT ripening process than in the standard ripening process, and one day and a half ealier in the TVA ripening process.

We can therefore conclude that the PT ripening process is very similar to the TVA one.

Conclusion

In this work, we used a parallel multi-objective evolutionary algorithm to model a cheese ripening process. Based on viability theory, the analysis yield a Pareto front of a set of viable trajectories. A analysis made by a cheese ripening expert allowed to select an interesting trajectory in this Pareto set. The major improvement of this optimal controlled trajectory is its shortening, which is in accordance with previous work on this topic. It has been experimentally verified that a 8 days trajectory (the TVA trajectory of figure 3 and4) was able to yield a similar cheese in terms of sensory panel [START_REF] Sicard | Development of a viability approach for reverse engineering in complex food processes: Application to a camembert cheese ripening process[END_REF] in comparison to the standard trajectory used in cheese ripening industry. Additionally, for the 8 days trajectory, the simulated process yield quantities that seems coherent with real experiments using expert-based optimised control settings. Further work on this topic will consist in using the Pareto optimal trajectory to control an experimental ripening process, in order to verify the precision and validity of the modeling.

Figure 1 :

 1 Figure 1: The viability graph calculated for 12 days of ripening. The time, respiration rate and cheese weigth viable value for each days are represented.

Figure 2 :

 2 Figure 2: Projection of Pareto front solutions using a principal component analysis (PCA).

Figure 3 :

 3 Figure 3: Comparison of cheese mass loss between PT trajectory (computed values), and TVA and SRT ones (measured values from real experiments in a ripening chamber).

Figure 4 :

 4 Figure 4: Comparison of respiration rates between PT trajectory (computed values), and TVA and SRT ones (measured values).

Table 1 :

 1 Parameters setting

	Population size

Table 2 :

 2 Summary of controls (temperature T • and relative humidity RH% ) applied to the optimised trajectory PT, the viable trajectory TVA , and a standard trajectory SRT.

tory (PT) for a 0.280 kg cheese. The controls of this trajectory are presented in table 2.