Huei-Fang Yang

Xavier Descombes

Charles Kervrann

Caroline Medioni

Florence Besse

Tracking Growing Axons by Particle Filtering in 3D + t Fluorescent Two-Photon Microscopy Images

Analyzing the behavior of axons in the developing nervous systems is essential for biologists to understand the biological mechanisms underlying how growing axons reach their target cells. The analysis of the motion patterns of growing axons requires detecting axonal tips and tracking their trajectories within complex and large data sets. When performed manually, the tracking task is arduous and time-consuming. To this end, we propose a tracking method, based on the particle filtering technique, to follow the traces of axonal tips that appear as small bright spots in the 3D + t fluorescent two-photon microscopy images exhibiting low signal-to-noise ratios (SNR) and complex background. The proposed tracking method uses multiple dynamic models in the proposal distribution to predict the positions of the growing axons. Furthermore, it incorporates object appearance, motion characteristics of the growing axons, and filament information in the computation of the observation model. The integration of these three sources prevents the tracker from being distracted by other objects that have appearances similar to the tracked objects, resulting in improved accuracy of recovered trajectories. The experimental results obtained from the microscopy images show that the proposed method can successfully estimate trajectories of growing axons, demonstrating its effectiveness even under the presence of noise and complex background.

Introduction

Analyzing how growing axons correctly reach their target neurons is essential for biologists to better understand the development of a nervous system. The advances in fluorescence imaging technology, such as two-photon microscopy [START_REF] Helmchen | Deep tissue two-photon microscopy[END_REF], have generated amounts of imaging data that enable the analysis of the dynamic behavior of developing neurons in real time and in three dimensions (3D). Analysis of the properties of axon growth requires detecting axonal tips that can appear as small bright spots and tracing their trajectories in low signal-tonoise ratios (SNR) image sequences. Figure 1 shows some examples of part of maximum intensity projection (MIP) images of the fluorescent two-photon microscopy volumes. As can be seen, in addition to high level of noise, the obtained volumes contain some bright spots that have appearances similar to the axonal tips and other biological structures that increase the difficulty in tracing the axonal tips. These factors make it difficult to manually identify the growing axonal tips and follow their trajectories. Besides, manual tracing by human experts to obtain the trajectories of axons is tedious and time-consuming. Therefore, developing automated tracking methods that produce reliable trajectories for further analysis is highly required.

Numerous methods have been proposed in the literature for tracking biological spot-like particles. A great majority of the tracking methods are based on two-step approach: the object detection/localization step and the data association step [START_REF] Genovesio | Multiple particle tracking in 3-D+t microscopy: method and application to the tracking of endocytosed quantum dots[END_REF][START_REF] Matov | Analysis of microtubule dynamic instability using a plus-end growth marker[END_REF][START_REF] Sbalzarini | Feature point tracking and trajectory analysis for video imaging in cell biology[END_REF]. The detection step aims at locating the objects of interest (e.g. the bright spots) by utilizing intensity thresholding [START_REF] Apgar | Multipleparticle tracking measurements of heterogeneities in solutions of actin filaments and actin bundles[END_REF], wavelet transform [START_REF] Olivo-Marin | Extraction of spots in biological images using multiscale products[END_REF], or Gaussian fitting methods [START_REF] Zhang | Gaussian approximations of fluorescence microscope point-spread function models[END_REF]. Since the detection step has been performed, the data association step solves the object correspondence problem between consecutive detection images [START_REF] Apgar | Multipleparticle tracking measurements of heterogeneities in solutions of actin filaments and actin bundles[END_REF]. Generally, these two steps are performed independently. Therefore, because the detection results may contain missing or incorrectly detected objects, the data association step tends to produce wrongly associated trajectories.

Recently, Bayesian tracking approaches have received much attention and have been successfully applied to spot tracking [START_REF] Cardinale | Bayesian image analysis with on-line confidence estimates and its application to microtubule tracking[END_REF][START_REF] Smal | Multiple object tracking in molecular bioimaging by rao-blackwellized marginal particle filtering[END_REF][START_REF] Smal | Particle filtering for multiple object tracking in dynamic fluorescence microscopy images: Application to microtubule growth analysis[END_REF]. One advantage of these approaches is that they integrate the detection and data association steps in a unified framework, which produces more accurate trajectories. Moreover, they can deal with nonlinear and non-Gaussian tracking problems. The Bayesian tracking is usually implemented by using a sequential Monte Carlo (SMC) method, that is the so-called particle filtering (PF) [START_REF] Isard | CONDENSATION -conditional density propagation for visual tracking[END_REF]. Particle filtering techniques use prior knowledge of the dynamics of the objects to predict the object state modes and use all the available information (noisy measurements) to compute optimally the posterior probability density function (pdf).

In this paper, we propose a tracking method, based on the particle filtering technique, to follow the trajectories of the individual growing axons in the fluorescent two-photon microscopy image sequences. The proposed method incorporates the important characteristics of growing axons. First, multiple dynamic models are used to better predict the locations of axonal tips. Furthermore, the proposed method incorporates information related to object appearance, motion characteristics of the growing axons, and filaments in the computation of the likelihood function. The integration of these three sources prevents the tracker from being distracted by other objects that have appearances similar to the tracked objects. It ensures that the tracker can successfully follow the traces of the individual axonal tips, resulting in obtaining more accurate trajectories.

The remainder of this paper is organized as follows. Section 2 presents an overview of the particle filtering framework, and Section 3 elaborates on the proposed particle filtering based tracking method that incorporates the characteristics of growing axons. The experimental results obtained from synthetic data and a 3D + t fluorescent two-photon microscopy data set are shown in Section 4, followed by conclusions and future work in Section 5.

Particle Filtering

The Bayesian tracking aims at recursively estimating the posterior probability density function, p (s t |z 1:t ), which describes the state s t at time t given a series of (noisy) measurements z 1:t from time 1 to time t. The estimation consists of two steps: prediction and updating.

Prediction The prediction step relies on the Chapman-Kolmogorov equation, that is,

p (s t |z 1:t-1 ) = p (s t |s t-1 ) p (s t-1 |z t-1 ) ds t-1 , (1)
where p (s t |s t-1 ) decides the prior dynamics.

Updating The update step uses the Bayes' rule to modify the prior density and obtain the posterior density when a measurement z t is available:

p (s t |z 1:t ) ∝ p (z t |s t ) p (s t |z 1:t-1 ) , (2)
where p (z t |s t ) represents the likelihood function (observation model). This recursive estimation of the posterior pdf is intractable, and the particle filtering is used as a numerical approximation. The particle filtering [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-gaussian Bayesian tracking[END_REF] uses a set of N s samples, each of which associates with a weight, s i t , w i t Ns i=1 , to represent the posterior p (s t |z 1:t ):

p (s t |z 1:t ) ≈ Ns i=1 w i t δ s t -s i t , (3)
where δ ( • ) is the Dirac delta function, and Ns i w i t = 1. These samples are propagated through time to approximate the posterior pdf at the subsequent steps via a proposal distribution q (s t |s t-1 , z t ). The importance weight of each sample, w i t , is updated as follows:

w i t ∝ w i t-1 p z t |s i t p s i t |s i t-1 q s i t |s i t-1 , z t . (4)
The final state at time t can be computed as the maximum a posteriori (MAP) estimator: ŝt = argmax

s i t p s i t |z t ∀i = 1 . . . N s . (5)
When applied to the tracking problems after a few iterations, the particle filtering has the degeneracy problem that only few particles carry significant weights. A resampling procedure is then required to generate a new set of equally weighted samples based on the importance weights.

Proposed Method

The section details the proposed method by starting with the description of the multiple dynamic models used in the proposal distribution. Then we present the observation model that integrates information related to object appearance, object motion, and filaments.

Multiple Dynamic Models

Before describing the dynamic models, the state of an axonal tip is first introduced. The state of an axonal tip at time t is represented by s t = (x t y t z t ẋt ẏt żt )

T , where x t , y t , and z t denote the position of the tip in the Cartesian coordinate and ẋt , ẏt , and żt represent the velocities in the x, y, and z directions. To better predict the position of the axonal tip at time t, the proposal distribution contains two dynamic models, defined as

q (s t |s t-1 , z t ) = π 1 p (s t |s t-1 ) + π 2 p (s t |s t-1 , s t-2 ) , (6)
where p (s t |s t-1 ) = N s t |s t-1 , σ 2 1 and p (s

t |s t-1 , s t-2 ) = N s t |s t-1 + v t-1 , σ 2 2
are the first-order Markov transition and the second-order Markov transition, respectively, where v t-1 is the velocity at time t-1, given by v t-1 = ( ẋt-1 ẏt-1 żt-1 ). The notation N :, σ 2 represents a Gaussian with variance σ 2 . The first-order Markov transition models the phenomenon that the growing axons sometimes stop growing for a short period of time and then resume the growth process, and the second-order transition captures the typical motion behavior of growing axons, which tends to follow the growing direction at time t -1. The weights π 1 and π 2 1 -π 1 control the mixture for each dynamic model.

Observation Model

The proposed observation model combines information about object appearance, object motion, and filaments in order to better follow the axonal tips.

Appearance Likelihood The appearance likelihood aims at modeling the prior knowledge about the image formation process in the microscopy imaging systems (point spread function (PSF)). Now, let I ideal t (x, y, z) be the ideal intensity, without the effect of noise, at point (x, y, z) at time t, and is given by [START_REF] Chenouard | Improving 3D tracking in microscopy by joint estimation of kinetic and image models[END_REF]

I ideal t (x, y, z) = i P i,t (x, y, z) + b t , (7)
where P i,t is the object profile, and b t is the background intensity. Here, Gaussian functions are used to approximate the PSFs of the microscopy. However, because of the noise generated during the image acquisition process, the observed intensity I obs t is the addition of the ideal intensity I ideal t and the noise N t that is assumed to be additive and Gaussian:

I obs t = I ideal t + N t . (8)
To compute the likelihood, the residual image between the observed image and the ideal image is first obtained by

R t = I obs t -I ideal t . (9)
After the residual image is obtained, the appearance likelihood is therefore computed as

p (R t |s t ) = 1 Z 1 exp (-α 1 ||R t || 2 ) , (10)
where α 1 is a constant, and Z 1 is a normalizing constant.

Motion Model It is established in biology that axons reach their target cells in the developing nervous system by the guidance of molecular gradients [START_REF] Goodhill | Predicting axonal response to molecular gradients with a computational model of filopodial dynamic[END_REF][START_REF] Krottje | A mathematical framework for modeling axon guidance[END_REF]. Accordingly, the growing axons do not abruptly change their growing directions.

In other words, the turning angles of the growing axons between consecutive time sequences are small, leading to smooth trajectories. The likelihood is designed to model this characteristic of the movement patterns of growing axons. Assume that r t is the motion direction vector at time t. The angle difference, d t , between the motion direction vectors at time t and at time t -1 is given as

D t = arccos (r t , r t-1 ) . (11)
Therefore the motion likelihood is defined as

p (D t |s t ) = 1 Z 2 exp (-α 2 D t ) , (12)
where α 2 is a constant, and Z 2 is a normalizing constant.

Filament Model Frangi's vessel enhancement filter [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] is applied to enhance the filaments in the images. The basic idea of Frangi's method is to use the eigenvalues of the Hessian (second-order information) to measure the vesselness. Now, let λ k represent the k-th smallest eigenvalue obtained from a Hessian matrix (i.e., |λ

1 | ≤ |λ 2 | ≤ |λ 3 |),
and the eigenvalues of an ideal tubular structure in a 3D volume has the following properties: (1)

|λ 1 | ≈ 0, (2) |λ 1 | |λ 2 |, and (3)
λ 2 ≈ λ 3 .
To quantitatively distinguish the tubular-like structure from blob-like and plate-like structures, two geometric ratios R A and R B are defined. They are

R A = |λ2| |λ3| and R B = |λ1| √ |λ2λ3|
, where R A is for distinguishing between platelike and line-like structures, and R B accounts for the deviation from a bloblike structure. Besides, because background pixels have small magnitudes of the eigenvalues, a measurement S is used to distinguish the tubular-like structures from background and is defined as S = j λ 2 j . Finally, a vesselness function V δ (p) that integrates the above 3 measurements, R A , R B , and S, at scale δ is given by

V δ (p) = 0 if λ 2 > 0 or λ 3 < 0 1 -exp -R A (p) 2 2γ 2 exp -R B (p) 2 2β 2 1 -exp -S(p) 2 2c 2 , (13)
where γ, β, and c are parameters that control the sensitivity of the line filter to the measures R A , R B , and S, respectively, and p denotes a voxel.

Equation ( 13) is analyzed at different scales to measure the vesselness in an image so the final estimate of the vesselness is the maximum filter response among all the scales, that is

V (p) = max δmin≤δ≤δmax V δ (p) , (14)
where δ min and δ max are the minimum and maximum scales at which relevant structures are expected to be found. and the others are suppressed when compared to the original MIP shown in Figure 2(a). Let l be the trajectory (a line segment) that an axon moves from time t -1 to time t and V l,t be the average value of the measured vesselness through this trajectory. The filament likelihood is therefore given by

p V l,t |s t = 1 Z 3 exp -α 3 1 -V l,t , (15)
where α 3 is a constant, and Z 3 is normalizing constant.

Joint Likelihood Assuming the independence of object appearance, object motion, and filament measurement, the joint likelihood for the observation model is defined as

p (z t |s t ) = p (R t |s t ) p (D t |s t ) p V l,t |s t . (16)

Experimental Results

The experiments were carried out on synthetic data and a fluorescent microscopy data set to evaluate the performance of the proposed method. The performance is quantitatively measured by comparing the computer generated tracking results to the ground truth based on evaluation metrics.

Experiments on Synthetic Data

Data To validate our method, a synthetic data set, with a size of 100 × 100 × 35, that simulates the movement of growing axons is generated. Different levels of Gaussian noise with standard deviation (STD) of 5 and 7 were added to the synthetic data, resulting in a total of 2 image sequences. Figures 3(a

Results

To quantitatively measure the performance of the proposed tracking method, the tracked trajectories were compared to manually generated ground truth by measuring the root mean square error (RMSE) [START_REF] Smal | Particle filtering for multiple object tracking in dynamic fluorescence microscopy images: Application to microtubule growth analysis[END_REF] of positions between the tracked trajectories and the ground truth. The RMSE is defined as

RMSE = 1 K K k=1 MSE 2 k , ( 17
)
where K is the number of independent runs, and

MSE 2 k = 1 M M m=1 1 T m t∈Tm r m,t -rk m,t 2 , ( 18
)
where r m,t is the true position of object m at time t, rk m,t is the estimated position, and T m is the number of frames where object m exists.

Figures 3(c) and 3(d) show the comparison between the tracking results generated by the proposed method, shown in red, and the ground truth, shown in green. By visual comparison, the computer generated tracks are consistent with the ground truth.

Using RMSE as the performance measure metric, the localization errors for the synthetic data set corrupted by Gaussian noise with STD = 5 and STD = 7 were 2.5565 and 2.6174 pixels, respectively, which indicates that the proposed method can correctly follow the movement of growing axons.

Experiments on Fluorescent Two-Photon Microscopy Data

Data A 3D+t fluorescent two-photon microscopy data set of the intact Drosophila fly brains was used to evaluate the performance of the proposed tracking method. The volume size at each time is 1012 × 548 × 16 (x × y × z) voxels, and 170 such volumes are generated. Each volume has a size of 133.75×72.36×12 µm 3 , with a resolution of 0.132 × 0.132 × 0.8 µm 3 /voxel. To analyze the dynamics of axons in the developing nervous system, the trajectories of the individual growing axons have to be tracked.

Preprocessing The fluorescent microscopy data set was obtained by imaging the living tissue at different time points so the subsequent volumes had to be aligned to obtain the correct motion patterns of growing axons. In addition, the raw images contain noise and complex background that have to be removed in order to obtain more accurate trajectories. Hence, prior to the tracking, two preprocessing steps, registration and background subtraction, were performed.

Registration Motion2D [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF], a software to estimate 2D parametric motion models in an image sequence, was used to register the images in order to eliminate the living sample's drift problem caused during the image acquisition process. To register the 3D+t image sequence, the MIP image at each time was first obtained, and the registration process was performed to estimate an affine motion model for each pair of consecutive MIP images. The obtained motion models were then used to register each individual frame. Though a 3D registration method that registers consecutive volumes may be more appropriate for this case, we observed that this 2D registration approach was able to give satisfying results.

Background Subtraction As demonstrated in Figure 1, the obtained images are noisy and contain other biological structures (background) that may affect the tracking performance. To remove noise and partial background, the HullkGround software [START_REF] Chessel | Computational geometry-based scale-space and modal image decomposition: application to light video-microscopy imaging[END_REF], which performs background subtraction by convex hull estimation, was used to decompose the registered data set into two dynamic components: a result, the image shown in Figure 4(b) still contains some parts of the background. Even so, the proposed tracking method still yielded better results on the moving objects images than the original images.

Tracking Results

The proposed tracking method was applied to the tracking of growing axons in the image sequence containing moving objects, which was obtained by background subtraction. The initial positions of the axonal tips were manually set by the user, and the tracking algorithm followed the traces of the growing axons in the subsequent images. In the experiments, a set of 200 samples was used in the particle filter. The weights (i.e. π 1 and π 2 in Equation ( 6)) for the two dynamic models were set to 0.5. The parameters α 1 , α 2 , and α 3 in the observation model were set to 0.2, 2.5, and 0.25, respectively. Figure 5 shows the 2D and 3D visual comparisons between the tracked trajectories produced by the tracking algorithm, shown in red, and the ground truth tracks, shown in green. As can be seen, the computer generated trajectories are consistent with the manually created ground truth in general, with minor differences between the estimated positions and the true positions, which is caused by the noise and complex background. This indicates that the proposed method was able to successfully follow the traces of axonal tips.

Using RMSE as the performance measure metric, for all the tracks, the overall localization error of the proposed method is 4.56 pixels.

The mean computation time using a Matlab implementation of the proposed method on a laptop with Intel Core i5 CPU (2.53 GHz) and 4GB memory was 1.8926 seconds for tracking one axonal tip at a time step. Note that this computation time does not include the time required for reading frames and applying Frangi's filter.

Performance Comparison Between Different Likelihood Combinations

To investigate how each likelihood term in the observation model affects the tracking results, we evaluated the tracking performance with 3 different likelihood combinations: (1) appearance and motion likelihoods and (2) appearance and filament likelihoods, and (3) appearance, motion, and filament likelihoods. The appearance likelihood is commonly used in most of the tracking methods and hence is retained in each combination. Hence, the main focus is on the investigation of how motion and filament likelihoods affect the performance. The localization errors for the 3 combinations in order were 7.51, 7.87, and 4.83 pixels. The integration of appearance, motion, and filament likelihoods in the observation model increases the localization accuracy.

Conclusions and Future Work

This paper presented a tracking method, based on the particle filtering technique, to follow the trajectories of the individual growing axons in fluorescent microscopy images. The proposed method uses multiple dynamic models in the proposal distribution in order to predict the positions of the axonal tips. Besides, information related to object appearance, motion, and filament is integrated in the computation of the likelihood. The fused information prevents the tracker from being distracted by the other objects that have appearances similar to the axonal tips; consequently, the proposed tracking method can successfully estimate the trajectories of the axonal tips in the image sequences that are noisy and have complex background. The experimental results obtained from a 3D + t fluorescent two-photon microscopy image sequence demonstrate the effectiveness of the proposed method.

Future work includes designing a detection method to automatically initialize the tracking process and quantitatively analyzing the tracked trajectories to gather statistical parameters of the growing axons. Current detection methods usually produce many false positives in detecting the axonal tips and thus a detection algorithm that is suitable for task is required. For analysis of the axons, the parameters, such as curvatures, turning angles, and growing lengths, will be extracted and used for mathematically modeling the growing patterns of axons. We expect that this statistical analysis and modeling will provide valuable insight into a better understanding of the characteristics of growing axons in the developing nervous system.

Fig. 1 .

 1 Fig. 1. Examples of part of maximum intensity projection (MIP) of the fluorescent two-photon microscopy volumes. (a) and (b) show part of MIPs from different times and demonstrate the dynamics of the growing axons. (c) and (d) are the close-ups of (a) and (b), respectively. The axonal tips are visualized as small bright spots, indicated by red circles. The 3D microscopy volumes are of low SNR ratios and contain other biological structures and some bright spots that have appearances similar to the axonal tips. These factors increase the difficulty in identifying and tracking the axonal tips.

Figure 2 (Fig. 2 .

 22 Fig. 2. MIPs of the original volume and of the Frangi's filtered volume. The tubular-like structures are enhanced and the other structures including background are suppressed in (b) after Frangi's filter is applied to (a). Note that the filter is applied to 3D volumes, and the MIPs of these 3D volumes are shown here for the visualization purpose.

  ) and3(b) show the generated synthetic images corrupted by Gaussian noise with STD = 5 and STD = 7, respectively.

Fig. 3 .

 3 Fig. 3. Generated synthetic images with different levels of Gaussian noise and tracking results. (a) and (b) show the synthetic images by adding Gaussian noise with STD = 5 and STD = 7, respectively. (c) and (d) show the tracking results in red and the ground truth in green. The results produced by the proposed method are consistent with the ground truth. Note that the contrast of images is enhanced for visualization.

Fig. 4 .

 4 Fig. 4. Obtained moving objects and partial background by using the HullkGround software. As can be seen, the original frame (a) contains noise and complex background, and the background subtraction step decomposes the original image into an image containing moving objects (b) and an image containing slightly moving background (c).

( 1 )

 1 an image sequence containing moving objects and (2) an image sequence showing the slightly moving background. Figures 4(b) and 4(c) show the image containing moving objects and the image containing the background, respectively. As can be seen, because of the high complexity of the background, only partial background is removed from the original image shown in Figure 4(a). As

Fig. 5 .

 5 Fig.5. Visual comparison between the tracking results of the proposed method and the manually created ground truth in 2D and 3D. The red trajectories are produced by the proposed method, and the green are the ground truth manually created by the expert; both are overlaid on the MIPs (a) and (b) and visualized in 3D (c). The computer generated tracks are consistent with the ground truth in general, with minor differences between the estimated positions and the ground truth positions. The differences are caused by the noise and by the effect of complex background.

Acknowledgement. This work was partially supported by ARC DATA (IN-RIA).


Jia Zhou
email: jia.zhou@cirad.fr

Christophe Proisy

Xavier Descombes

Josiane Zerubia

Guerric Le Maire

Yann Nouvellon

Pierre Couteron

Tree crown detection in high resolution optical images during the early growth stages

Keywords: tree detection, Eucalyptus plantation, marked point process, multi-date detection I

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Tropical plantation forests provide an important share of the global wood supply, where Eucalyptus is the most widely planted hardwood genus. In Brazil, eucalyptus plantations occupied a large surface of about 4.3 million ha in 2008. These fast-growing trees have a rotation length of typically 5 to 7 years [START_REF] Le Maire | Leaf area indix estimation with MODIS reflectance time series and model inversion during full rotations of Eucalyptus plantations[END_REF]. Their growth undergoes two periods: expansion of crowns and heights in the first year and development of tree heights in the following years. Monitoring work of this study is mainly realized in the first year.

Monitoring these plantations on large surfaces is rather difficult and requires a lot of trained manpower. The canopy high-resolution images provided by modern satellite-borne remote sensing techniques give another representation of extensive plantations in which individual trees may be documented via the monitoring of the tree crowns. The full uses of these images would help to efficiently analyze and supervise the planting work. The monitoring of the plantation during the first year raises specific questions concerning the detection of trees that die or for which growth is insufficiently fast, because the corresponding locations are to be replanted. It is also important to map areas of bad growth for elucidating the causes of failure. Documenting how the size of the plantation that grows either faster or slower than the average is an indirect mean to define the fertility of the land and the relative performance of the clones and strains that are used.

To do so, we employ the marked point process [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF][START_REF] Zhou | Tree crown detection in high resolution optical and LiDAR images of tropical forest[END_REF] methods to identify every eucalyptus crown.

II. FIELD INVENTORY

The Eucalyptus in this study were planted in November 2009, with a density of about 1670 trees/ha, in São Paulo State, south-eastern Brazil. Over a large area (>200 ha), they were planted in rows, with a distance of about 3 meters between two lines and a distance of about 2m (± 0.5m) between two trees in a line. We selected 4 plantation plots, where field measurements were made concomitantly with the acquisitions of satellite images: in May and August 2010. In each plot, 84 trees were chosen (14 consecutive trees in a line over 6 consecutive lines) for measurement of the crown diameters in line and inter-line directions, and the distance between two adjacent trees in a line. These data provided ground-truth to assess the crown detection performance. The crown centers of these trees were visually localized in the images by an expert taking into account the external information of approximate tree density and spacing.

III. VERY HIGH SPATIAL RESOLUTION IMAGES

The two images provided by the WorldView-2 (WV2) satellite are panchromatic and of the spatial resolution: 0.5m. One was acquired on 11th May 2010, centered at (23.031° S, 48.694° W), with the off-nadir angle: 19°, and the target azimuth: 233°. The second was acquired on 1st August 2010, centered at (22.997° S, 48.685° W), with the off-nadir angle: 7°, and the target azimuth: 114°.

Extracts of those images covering the same zone of the monitored plantation at the two successive dates are shown in Fig. 1:

A. Marked point process

We consider a multi-date marked point process to model the plantation both in May and in August [START_REF] Zhou | Tree crown detection in high resolution optical and LiDAR images of tropical forest[END_REF]. An object in the configuration is defined by a pair of disks (c, r (1) , r (2) ), where c є P is the center, r (1) is the radius of the tree crown in May and r (2) in August. We assume that r (2) ≥r (1) . The object space is given by: where X M and Y M are the width and length of the images I 1 and I 2 and: r m (respectively r M ) is the minimum (respectively maximum) radius of the tree crown.

The marked point process can be expressed by an energy decomposed as the sum of a prior and a data term: where is a configuration.

We consider a prior penalizing overlaps between neighboring objects: where is the number of points of intersecting disks in the image I 2 (note that the disks are bigger in the second image).

To define the data energy term, we consider an object based model as follows:

For each image I, we compute a statistical test Test j (c i , r i (j) ).To compute the function Test j (c i , r i (j) ), we consider the set of pixels inside the disk and compute their mean and variance . We also compute the mean and variance of pixels inside the concentric crown of one pixel width.

If , then Test j (c i , r i (j) ) = 1, which means that the disk does not fit a tree, assuming that the tree radiometry is higher on the ground. Otherwise, we compute a Student test:

We then compute the p-value p j (c i ,r i (j) ) associated with this statistical test and consider a threshold t, as a parameter. Then the data term is written as follows:

Therefore, an object fits the data if the p-value is smaller than 1-t.

To use the information of both images, we assume that a tree should be detected if the test function is negative in at least one image. Therefore, we define the data term as follows:

To optimize the model, we consider a simulated annealing based on a multiple birth and death process (MBD). This process has first been proposed in [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF].

The algorithm simulating the process is defined as follows:

1. Initialize the temperature parameter T = 1, the birth rate

and alternate birth and death steps.

2. Birth step: for each pixel s in the image I, if there is no disk centered in s, we add an object in s, with a random radius between r m and r M , with probability B(s), where B(s) is proportional to δ and the data term in s (see [START_REF] Descamps | Flamingo detection using marked point processes for estimating the size of populations[END_REF] for details).

3. Sorting step: once the birth step is finished, we compute the data term Test(w i ) for each object in the current configuration of objects. Then, we sort them, in decreasing order, according to their data energy.

4. Death step: for each object xi taken in this order, we compute the death rate as follows:

, where :

and being the number of disks in W intersecting w i .

5. Convergence: if the process has not converged, which means than the configuration has changed during the last iterations, then decrease T and δ, according to a law defined in [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF], and go back to step 2.

B. Detection score

To assess the detection results in tree numbers and positions compared to the field reference, we used a score which was already introduced in [START_REF] Larsen | Comparison of six individual tree crown detection algorithms evaluated under varying forest conditions[END_REF] for comparing different tree crown detection algorithms. For this, a matching work has to be made first between the detected disk set and the reference tree set, then the successful detection number, the omission and false detection number can be calculated. For a reference tree be considered to be detected by a disk (a successful detection), the disk center should fall within the extent of the corresponding tree crown. A reference tree can be matched with only one disk, and other disks not matched are considered as false detections. A tree not matched with any disk is considered as an omission. The detection score is defined as follows:

where Ns is the successful detection number, No is the omission number and Nc the false detection number.

C. Interpolation

The sizes of planted Eucalyptus crowns are first limited by the space between adjacent trees in the same line because of competition for spaces. The limitation between two lines occurs later. The diameters in the inter-line direction begin to exceed that in line direction in the first year of plantation. Nevertheless, the mean diameters (the arithmetic means of diameters in two directions) exceed rarely 2.7m. Furthermore, the field measurements of diameters were based on the branches which extend the most from the crowns center in the line and inter-line directions, and may often be not visible in the images. Thus, in the 0.5m resolution images, a tree crown covers just a few pixels, with a radius length of 1 or 2 pixels. These small crowns make it difficult for the marked point processes to give correct detections both in position and in crown size. We were working in the limited conditions for these young populations. Better resolution images are needed for more detailed analysis.

A solution to this problem is the resampling of the original images by interpolation. The usual interpolation methods [START_REF] Thévenaz | Image Interpolation and Resampling[END_REF] for images like nearest-neighbors, bilinear, and bicubic were tested at different scales between 2 and 10. The detection works better on the interpolated images, usually when they were resampled at 5, which means at the resolution of 0.1m. The detection results were little influenced by the choice of interpolation method. Therefore, we will just show the results on the bilinear interpolated images in the following.

V. RESULTS

The detection results were first analyzed in terms of tree number, Table 1 contains the list of statistics in each plot. On the whole, the detection is encouraging for tree number counting. For each plot, almost all the reference trees were detected (between 89% and 95%) with just a few false detections (between 0 and 1.3%). The detection in crown position and delineation of one plot are shown in Fig. 2. For the crown size estimation, we compared, for every detected tree in the reference zones, the estimated average crown diameter given by the process with the measured crown diameter in the field. In Fig. 3, we can see these comparisons of diameters in 4 plots. Note that the precision of estimated crown diameters depends on the image resolution, so the minimum difference between successive diameter values is twice the spatial resolution of the image interpolated with bilinear method. However, without taking this into account, many points in the scatter graphs were still away from the line x=y. The dashed lines built an interval of [-0.5m, 0.5m] around the line x=y, with 0.5m corresponding the image resolution.

When considering the scatter graph of August, the clouds of points were moved a little to the right, which means that in these plots, the detection process tended to underestimate the crown sizes. This phenomenon is related, on one hand, to the fact that the measurement of diameters in the field was based on the branches which extended the most from the crown centers, but were not always detectable on images, as we mentioned earlier; and on the other hand, to a weakness of the detection model, in which a tree crown is simply modeled by a disk. Indeed, real tree crowns are not perfect disks, thus crowns touched each other along the line direction and extended more in the inter-line direction, which was the case in August. The penalization of large overlapping of disks in the detection algorithm, with purpose of avoiding false detections, disfavors big disks in this case. The diameters of big disks are usually limited by the distances between adjacent trees in a line (if the crowns are touching each other in the line direction), which were about 1.9m according to field measurements, while the mean of measured diameters were about 2.1m in August. From the detected crown diameters in May and in August, the crown growths were deducted for individual trees, which are represented in Fig. 4 as histograms. Compared with the field measurements in Fig. 5 (with the modes at 0.4~0.5m), we note that these growths were widely underestimated (with the modes at 0.2~0.3m in Fig. 4). However, it was impossible to get a better accuracy on crown diameter estimation with these two images which are limited in spatial resolution, so the estimations of crown growths were also influenced. In this study, we were in the limited case of the young Eucalyptus stage analysis (owing to the ratio of crown sizes and the image resolution), using satellite images of spatial resolution 0.5m. Although interpolation methods could improve detection results, they were not the solution to the essential problem of resolution limit, because noise might be introduced during the interpolation process. The detected crown sizes were of the same order of magnitude as the measured sizes in the field, however, an accuracy of better than 0.5m cannot be asked. During the period between the two images in this study (81 days), the best crown growth was about 0.7m in diameter, which was not detectable with these images. But, the detection results remain consistent with the dynamics of the young eucalyptus population, which shows the potential of our detection method.

Figure 1 .

 1 Figure 1. Image extracts of the same zone at two dates (May, August)

Figure 2 .

 2 Figure 2. Detection result in May (on the top) and in August (at the bottom) of one plot. The points represent detected tree positions, and the circles represent delineated crowns.

Figure 3 .

 3 Figure 3. Comparisons of the average diameters of the detected tree crowns with field measured crown diameters in May (on the top) and in August (at the bottom) of 4 plots.

Figure 4 .

 4 Figure 4. Distributions of detected growths in 4 plots.

Figure 5 .

 5 Figure 5. Distributions of measured growths in 4 plots.

TABLE 1 .

 1 Statistics of detection results in 4 plots

	Plot	Tree number	Detected tree number	Successful detections	% detected trees	False detections	% false detections	Detection scores
	1	83	75	74	89.2	1	1.3	88.1
	2	84	79	79	94.0	0	0	94.0
	3	84	77	77	91.7	0	0	91.7
	4	83	80	79	95.2	1	1.3	94.0
		VI. DISCUSSION

We also tried to apply the method on other population of tropical forest stands, as natural mangrove forest in French Guiana, where detections were much more difficult owing to the complexity and high vegetation density of the population. In addition, to better test and adjust the new method, we used the simulated optical DART (Discrete Anisotropic Radiative Transfer) images [11,12], where exact field inventory could be provided for detection analysis.

ACKNOWLEDGMENT

These satellite images were acquired within the ORFEO program, a CNES Accompaniment Program for the use of PLEIADES images. We acknowledge this project for funding the images, and Claire Tinel (CNES) for her help with the images acquisitions. We are grateful to the Eucflux project (http://www.ipef.br/eucflux/) for providing the inventories of the eucalyptus plantation plots.