Claire Coiffard-Marre

Ségolen Geffray
email: segolen.geffray@math.unistra.fr

Object detection and automatic measurements with a data-driven algorithm for marked point process optimization Long technical extended version

Keywords: Spatial marked point processes, maximum a posteriori, reversible jump Markov Chain Monte-Carlo, Simulated annealing, Radon transform, circular Hough transform

This paper deals with object detection and automatic measurements from optical microscopy and scanning electron microscopy (SEM) 2-dimensional images. More specifically, we consider the problem of detecting and measuring circular cells and rectilinear thick fibers. We adopt the spatial point process viewpoint and follow the maximum a posteriori (MAP) principle to obtain an optimal object configuration. We design a novel data-driven version of the data term in the MAP criterion in the fiber case, of the reference spatial marked point process and of the Markov Chain Monte-Carlo (MCMC) sampler, both in the cell case and in the fiber case. To this purpose, information from the image under consideration is incorporated by means of either the Radon transform or the circular Hough transform. We apply the proposed method to a real cell optical microscopy image and to a real fiber SEM image. Our results show that the proposed algorithm is fast, reliable and stable. We also produce the required characterization of the object morphology.

Introduction

Images are a major source of information when dealing with microscopic objects. For example, think of cells or of polymeric fibers which are roughly of the order of the micrometer. In polymer engineering, systematic measurements are made repeatedly from polymer scanning electron microscopy 1 (SEM) images, as outlined in Ahirwal et al. [START_REF] Ahirwal | From Self-assembly of Electrospun NanoFibers to 3D cm-thick Hierarchical Foams[END_REF], Lavielle et al. [START_REF] Lavielle | Controlled formation of polycaprolactone ultrathin electrospun nanofibers by hydrolytic degradation approach[END_REF] or Shin et al. [START_REF] Shin | Determination of electrospun fiber diameter distributions using image analysis processing[END_REF] for morphology evaluation purposes. In cell biology, the authors of Vigneron et al. [START_REF] Vigneron | Adaptive filtering and hypothesis testing: Application to cancerous cells detection[END_REF] developed an automated low-level method of cell detection and counting from optical microscopic images. With these two examples in mind, a crucial goal for biomaterial engineering consists of automatically detecting such objects in a given image and of being able to characterize their location and morphology. Object detection from image data can be done using different techniques. Low-level image analysis methods present the major drawback of being sensitive to local image imperfections such as noise, poor illumination, blur... Moreover, once the detection is obtained with such a low-level image analysis method, it may be tedious to deduce object measurements. As a consequence, we turn to high-level methods in which an image I is considered as an imperfect (noisy, blurred, non-uniformly illuminated,...) observation of a random process X in an observation window X . When the objects of interest can be modelled as simple parametric shape such as circles, squares, triangles, thick straight lines, then a suitable class of random processes for that modelling purpose is the class of spatial marked point processes. We refer the reader to the books by Møller and Waagepetersen [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF] or by Van Lieshout [START_REF] Van Lieshout | Markov point processes and their applications[END_REF] for a detailed presentation of spatial point processes, their sampling and practical applications. Object detection and counting can be made by fitting a marked point process to the image as largely demonstrated in the literature. We recommend the book by Descombes et al. [START_REF] Descombes | Stochastic geometry for image analysis[END_REF] and the numerous references therein for applications of spatial point processes to practical image analysis situations. This paper is concerned with providing an accurate, reliable and robust determination of cell and fiber location and morphology from microscopy images within a reasonable computation time and complexity. Specifically, we aim at detecting cell and fiber objects from microscopy images and at estimating the cell number and the cell radius distribution on the one hand, and the fiber width and orientation distributions on the other hand. We model cells as circles and fibers as thick straight lines within the framework of an appropriate spatial marked point process. We adapt the marked point process method presented by Descombes et al. [START_REF] Descombes | Stochastic geometry for image analysis[END_REF] for its ability to provide straightforward automatic measurements once the object detection is done. Two kinds of inferential approaches may be distinguished when fitting a spatial marked point process to a given image, the Bayesian approach on the one hand, and the so-called detector approach according to by Descombes et al. [START_REF] Descombes | Stochastic geometry for image analysis[END_REF] on the other hand. The Bayesian approach has been used successfully in a number of practical situations and recently in Konomi et al. [START_REF] Konomi | Bayesian object classification of gold nanoparticles[END_REF]. However, the Bayesian method requires to specify a likelihood for the grey level intensity of pixels within and outside the objects. This makes the Bayesian method very sensitive to local defects in the image such as noise, blur or illumination gradient and also to the fact that grey level intensities may vary from on object to another which is the case in practice (see Figure 1). For these reasons, we adopt the detector approach as inferential method. Instead of specifying a likelihood for the grey level intensity of pixels within and outside the objects, one only specifies a so-called data term which accounts for the plausibility of a given object configuration with respect to the image at hand. The detection then consists of finding the optimal object configuration which is done by maximizing an objective function. Due to the intractability of the resulting nonconvex optimization problem, a simulated annealing algorithm (Kirkpatrick et al. [START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF]) together with a Reversible Jump Markov Chain Monte-Carlo (RJMCMC) sampler (Green [START_REF] Green | Reversible jump Markov Chain Monte Carlo computation and Bayesian model determination[END_REF]) are used. This method has demonstrated its wide applicability. However, the algorithm convergence to equilibrium may be really slow and the convergence may be somewhat difficult to diagnose. To speed up the algorithm, Brooks et al. [START_REF] Brooks | Efficient construction of reversible jump MCMC proposal distributions (with discussion)[END_REF] developed a general method for finding efficient proposals for RJMCMC algorithms. Their proposals are constructed from Taylor expansions of acceptance probabilities with the drawback of being computationally demanding. More specific data-driven versions for the RJMCMC algorithm have already been considered by Tu and Zhu [START_REF] Tu | Image segmentation by data-driven Markov Chain Monte Carlo[END_REF]. Their paper is segmentation oriented -and also Bayesian oriented. Their method transposed to our cell and fiber detection context would use as proposal distribution a finite mixture of Gaussian distributions the components of which are centered at the local modes of the Radon or Hough transform. Though interesting, their method increases the number of unknown model parameters since numerous additional variance components have to be calibrated making the calibration step tedious.

We present a novel data-driven version of both the objective function in section 2 and of the optimization algorithm in section 3. To this aim, information from the image is incorporated by means of either the Radon or the circular Hough transform. We show in section 4 that our method, though quite simple, leads to improved convergence properties of the algorithm while producing accurate, reliable and robust results. Some technical elements proving the convergence of the proposed algorithm to the target are postponed in the Appendix.

Objective function to be optimized

Object detection is obtained by computing the optimal object configuration according to the MAP criterion:

x m = argmax x m ∈Ξ P (I|X M = x m )P (X M = x m ). (1)
The first term P (I|X M = x m ) is the so-called data term and represents the probability of obtaining image I under the spatial marked point process model specified for the spatial marked point process denoted as X M and detailed below. The second term P (X M = x m ) is the so-called a priori term and represents the probability of obtaining the object configuration x m as a realisation of the random process X M . Notice that the target distribution of which we are looking for global maxima is the a posteriori distribution hence the denomination of Maximum A Posteriori principle.

We first design the a priori term as a spatial marked point process used to model cells or fibers in subsection 2.1. Then, in subsection 2.2, we design the data term according to the detector approach. We then derive the objective function according to the maximum a posteriori principle. Our main preoccupation is to incorporate information from the image data to make the algorithm efficient, robust and reliable.

A priori term

We first introduce some useful mathematical formalism. Let λ Leb be the Lebesgue measure in R 2 and let B be the associated Borel σ-algebra on R 2 . We assume that the location space X is a compact subset of R 2 with positive measure. Let X be a spatial point process on X so that on can think of X as a random countable subset of X . A point configuration x is an unordered sequence of locations in X . Let us denote as n(x) the cardinality of x. When n(x) = n for some n in N * , then the configuration x may be written as x = {x 1 , ..., x n }. We set (X , B X , λ Leb X ) to be the natural restriction to X of (R 2 , B, λ Leb ). In practice, the location space X is given by the image boundaries that is X = [0, x max ] × [0, y max ].

We then assume that the mark space M is a compact subset with positive measure of R d for some d ≥ 1. The mark space M is supposed equipped with the corresponding σ-algebra F M and a probability measure P M having density p M with respect to λ Leb M the restriction of Lebesgue measure on R d to M.

We introduce a marked point process on X × M as

X M = {X M i i := (X i , M i )} i=1,...,n : {X i } i=1,...,n ∈ X, M i ∈ M, i = 1, ..., n
meaning that a random mark M i is attached to each random point X i . A marked point is also termed as an object.

The mark space M naturally depends on the objects geometry. Cells are modelled as circles while fibers are modelled here as thick straight lines. In the cell case, the parametrization is straightforward. The object location is given by the circle center x = (u, v) in [0, x max ] × [0, y max ] and the associated circle radius is given by r in [0, r max ]. In the fiber case, the object location is given by (α, p, ε) in X as defined on figure 2, and the mark of the fiber is w in M = [w min , w max ] which corresponds to the fiber width.

For n in N, let Ξ n be the set of all unordered configurations x m = (x, m) = {(x 1 , m 1 ), ..., (x n , m n )} consisting of n marked points (x i , m i ) in X ×M with Ξ 0 being the empty configuration. For simplicity, we denote as n(x) the cardinality of a configuration x m since the number of marked points corresponds to the number of unmarked points. The configuration space can be written as Ξ = ∞ n=0 Ξ n and is equipped with the σ-algebra F Ξ defined by the mappings

ϕ A,n : {x m 1 1 , ..., x mn n } → n i=1 I(x m i i ∈ A)
that count the number of marked points in Borel sets A ⊆ X × M.

α, ε p w • α ∈ [0, π] • p ∈ [0, d max (α)] • ε ∈ {-1, 1} • w ∈ [w min , w max ]
Figure 2: Fiber parametrization

The Poisson process is a model for "no interaction" situations or of complete spatial randomness since for any disjoint subsets A and B of X , the events {X ∈ A} and {X ∈ B} are independent. In many applications, the reference Poisson process is the homogeneous Poisson process with independent marks having uniform distribution. It turns out that it is possible to get information about the likely localization and geometry of the objects in an image by using some image processing techniques for example. It seems relevant to use this information within the reference process definition, as in AlAwadhi et al. [START_REF] Al-Awadhi | Statistical image analysis for a confocal microscopy two-dimensional section of cartilage growth[END_REF]. We follow the approach of Baddeley and Van Lieshout [START_REF] Baddeley | Stochastic geometry model in high-level vision[END_REF] and Van Lieshout [START_REF] Van Lieshout | Markov point processes and their applications in high-level imaging[END_REF] who first define a reference measure denoted as here as π Λ that incorporates some features of the object distribution but not their mutual interactions and then set the full model as Gibbs point process defined by a Radon-Nikodym derivative with respect to π Λ . The point process taken as reference process hereafter is the inhomogeneous marked Poisson process defined for A ∈ F Ξ by

π Λ (x m ) = n∈N exp(-Λ(X )) n! (X ×M) n I({(x 1 , m 1 ), ..., (x n , m n )} ∈ A) n i=1 λ(x i )p M (m i )dλ Leb X (x i )dλ Leb M (m i ) (2)
where the integral term for n = 0 reads as I(∅ ∈ A), where Λ is the intensity measure, where

λ(•) = dΛ(•)
dλ Leb X is the intensity function and where p M is the density of the distribution P M on M with respect to λ Leb M . This says that the measure π Λ distributes the object locations independently in X according to a Poisson distribution with intensity measure Λ(.) and that the corresponding shapes are chosen in M independently according to the distribution P M . We design the intensity measure Λ and the mark distribution P M so that they incorporate some knowledge from the image data according to some deterministic image analysis methods. We set the intensity measure Λ to be equal to the image circular Hough transform in the cell case and to the image Radon transform in the fiber case.

In the cell case, we set the mark distribution P M to be equal to the radius distribution calibrated from the image in the following way. We compute the image circular Hough transforms for radii within a prespecified range [r min , r max ]. We then extract the local maxima which correspond to potential circles. We calibrate the radius distribution as the distribution of the radii associated with the obtained potential circles.

We planned to calibrate the fiber width distribution in a similar way. However, in practice, the number of fibers detected from our image by image analysis techniques is too small to derive a corresponding width distribution in a reliable way. Consequently, we kept the fiber width distribution as uniform within a prespecified range [w min , w max ].

To emphasize the difference with the usual choice, we note that the unit rate homogeneous marked Poisson process with uniformly distributed independent marks is usually chosen as reference process which boils down to taking λ(.) ≡ 1 and p M (.) = I(. ∈ M) in equation [START_REF] Al-Awadhi | Statistical image analysis for a confocal microscopy two-dimensional section of cartilage growth[END_REF]. Another justification for our reference process choice will also appear in the next section.

Interactions between objects have to be taken into account in order to prevent from excessive object superposition. Contrarily to Poisson processes, Gibbs point process allow the modelling of dependence of spatial point patterns. They can be constructed by specifying a Radon-Nikodym derivative with respect to π Λ . A repulsive pairwise interaction model which penalises overlap between objects is appropriate for our applications. These processes have densities proportional to

p(x m ) = x m i i ∈x m β(x m i i ) x m i i ,x m j j ∈x m :x m i i =x m j j ϕ(x m i i , x m j j ) with β : X → [0, ∞) and ϕ : X ×X → [0, ∞). If the interaction function ϕ satisfies ϕ ≤ 1, the process is repulsive and if in addition X ×M β(x m 1 1 )λ(x 1 )λ Leb X (dx 1 )p M (m 1 )λ Leb M (dm 1 )
< ∞, then the process is locally stable which implies that the process is well defined and satisfies the hereditary property

(that is p(x m ) > 0 =⇒ p(x ′ m ′ ) > 0 , ∀x ′ m ′ ⊂ x m ).
The local stability condition will also be of use to obtain the convergence properties of the Monte-Carlo dynamics in section 3. For both cell and fiber images, we choose an area-interaction process. This choice is quite classical in the literature, see eg Descombes et al. [START_REF] Descombes | Stochastic geometry for image analysis[END_REF] and the references therein. This process corresponds to the choices β(.) ≡ β for a β > 0 and ϕ(x

m i i , x m j j ) = exp -γ p t x m i i , x m j j
with γ p > 0 and with t x m i i , x m j j representing the overlap degree between the two objects x m i i and x m j j in the following way

t x m i i , x m j j = Area x m i i ∩ x m j j min Area (x m i i ) , Area x m j j
.

The area interaction process then has the following unnormalized density with respect to π Λ :

h p (x m ) = β n(x m ) e -γps(x m ) , (3)
where

s(x m ) = i<j t x m i i , x m j j
stands for the overall overlap degree in the current object configuration.

Data term determination

We now deal with the so-called data term. According to the detector approach, its expression is given by

h d (x m ) = exp   -γ d x m i i ∈x m v d (x m i i )   (4)
where γ d > 0 is a weight parameter and where v d : X × M → R is a so-called potential function which quantifies the confidence one may have in the fact that x m i i is an object actually present in the image. With such a specification, the plausibility of a configuration is sort of quantified by the addition of the plausibility of each object in the current configuration. The potential function in x m i i is classically given by some statistical distance between the gray levels of the pixels within the object x m i i and of those in its neighbourhood. The statement leading to this definition is twofold. Firstly gray levels are homogeneous inside an object. Secondly gray levels inside an object contrast with grey levels in the object neighbourhood. To assess the distance between the grey levels inside and in the neighbourhood of the object, generally either the Bhattacharaya distance is used (Perrin et al. [START_REF] Perrin | A Non-Bayesian Model for Tree Crown Extraction using Marked Point Processes[END_REF]), or a test statistic such as Student's test statistic (Lacoste [START_REF] Lacoste | Road network extraction in remote sensing by a Markov Object process[END_REF]) or Kolmogorov-Smirnov's test statistic (Imberty and Descombes [START_REF] Imberty | Simulation de processus objets : Etude de faisabilité pour une application à la segmentation d'image[END_REF]). The potential function then aims at penalizing objects for which the preceding distance is too small. We consider here separately the fiber case and the cell case. In the cell case, we adopt the potential function used in Perrin et al. [START_REF] Perrin | A Non-Bayesian Model for Tree Crown Extraction using Marked Point Processes[END_REF]. Let us denote as x m i i the circle in configuration x m under consideration and let b i (ρ) be its neighbourhood defined as the outer pixels no distant than more than ρ pixels from the object boundary. As explained above, the statistical distance between the grey levels inside object x m i i and in its neighbourhood is set as the Bhattacharyya distance:

d b (x m i i , b i (ρ)) = ( m 0 -m b ) 2 4 σ 2 0 + σ 2 b - 1 2 log   2 σ 2 0 σ 2 b σ 2 0 + σ 2 b   (5)
where m 0 and m b are the means of the grey levels respectively inside the object x m i i and in its neighbourhood b i (ρ), where σ 2 0 and σ 2 b are the variances of the grey levels respectively inside the object x m i i and in its neighbourhood b i (ρ). The potential function is then defined as

v d (x m i i ) =    1 - d B (x m i i ,b i (ρ)) d 0 1/3 if d B (x m i i , b i (ρ)), exp - d B (x m i i ,b i (ρ))-d 0 d 0 -1 otherwise. ( 6
)
In the fiber case, we propose a new expression of the potential function which relies on the Radon transform. Indeed, up to our knowledge, there does not exist any implementation of spatial marked point process for thick straight lines detection. One possibility would consist of using a potential function expression similar to that used in the cell case or to that used in the segment detection case by Lacoste et al. [START_REF] Lacoste | A Comparative Study of Point Processes for Line Network Extraction in Remote Sensing[END_REF]. A major drawback is the exploding computation time, since the thick straight lines have a very large number of pixels. Moreover, the presence of intersecting fibers in the neighbourhood of the fiber under consideration entails that the grey level intensities of the pixels in the neighbourhood of the considered fiber are inhomogeneous which makes irrelevant the use of (5) or of any similar distance between the grey levels inside and outside the fiber. Let us denote as x m i i the object in configuration x m under consideration. Let RT (x i ) be the value of the image Radon transform at the point having location x i . We recall that the Radon transform is defined as the integral of the image grey levels intensity over straight lines. Consequently, the larger RT (x i ) is, the more likely there is a line at location x i . To incorporate some additional information about the mark m i and get a more precise potential, we design a second term that accounts for the fact that the fibers are in light grey on a dark grey background. This term is defined as

k(x m i i ) = 1 n p (x i ) max (1, #{pixels of thick line x m i i having grey level intensity < τ }) , (7)
where n p (x m i i ) is the number of pixels in x m i i and where the threshold τ needs to be calibrated. Then, the potential v d (x m i i ) is obtained through a linear transform of

h i = RT (x i ) k(x m i i )
involving three additional parameters s, s ′ and V max given by

v d (x m i i ) =      V max if h i < s, 1 -2 h i -s s ′ -s if s ≤ h i ≤ s ′ , -1 otherwise. ( 8
)
As a conclusion to this section, we are in position to derive the objective function to be optimized. The data term P (I = y|X M = x m ) is taken proportional to the unnormalized density h d (x m ) being defined by equation ( 4). The prior term P (X M = x m ) is taken proportional to the unnormalized density h p (x m ) being defined by equations ( 3) and [START_REF] Imberty | Simulation de processus objets : Etude de faisabilité pour une application à la segmentation d'image[END_REF]. With these statements in mind, the optimization problem in (1) may be rewritten as

x m = argmax x m h(x m ), (9)
where h(x m ) = h d (x m )h p (x m ) is the objective function according to the MAP principle. The next section explains how to find a configuration maximizing the unnormalized density h.

Optimization

The task corresponding to the non-convex optimization problem in ( 9) consists of locating a global mode of h(.). For that purpose, we use the simulated annealing algorithm. The simulated annealing (SA) algorithm produces a random sample from h(.) which is biased overwhelmingly in favour of the most probable value x m . Specifically, the SA algorithm consists of repeatedly sampling configurations under the process having unnormalised density h 1/T (.) with T going to 0 as the number of iterations grows. Each distribution proportional to h 1/T (.) has its mode at x m , and as T decreases the mode becomes more and more exaggerated. Consequently, if we generate a random sample from each successive distribution, there is an increasing probability of producing the target x m . We adopt the classical geometrical cooling scheme which boils down to considering a positive sequence of temperatures (T k ) k∈N with T k = cT 0 where T 0 stands for the initial temperature and where 0 < c < 1. Both T 0 and c need to be calibrated. The simulated annealing algorithm is presented in algorithm 1. It involves a sampling step. Gibbs point processes are generally intractable because of the unknown normalizing constant of the likelihood. Therefore the sampling step of the SA algorithm motivates the use of a reversible jump Markov chain Monte-Carlo (RJMCMC) sampler also termed as Metropolis-Hastings-Green sampler. The construction of this sampler was given by Green [START_REF] Green | Reversible jump Markov Chain Monte Carlo computation and Bayesian model determination[END_REF] as an extension of the Metropolis-Hastings algorithm to general state-space. We refer the reader to the monographs by Robert and Casella [START_REF] Robert | Monte Carlo Statistical Methods[END_REF] for a detailed presentation of MCMC methods. Run the sampler with unnormalised density h(.) 1/T

4:

Change T ← cT with c < 1 Consequently, the optimization problem (9) moves to building an algorithm that simulates a Markov chain having as equilibrium distribution the measure π of our spatial marked point process defined as

π(A) = A h(x m )π Λ (dx m ), A ∈ F Ξ
and the realizations of which are configurations of marked points of X × M. We build the Markov chain to be π-invariant, aperiodic and geometrically ergodic, as proved in the appendix.

The general Metropolis-Hastings framework simulates Markov chains induced by the transition kernel

P (x m , A) = A α(x m , x ′ m ′ )Q(x m , dx ′ m ′ ), x ∈ Ξ, A ∈ F Ξ , x / ∈ A and P (x m , {x m }) = x ′m ′ =x m (1 -α(x m , x ′ m ′ ))Q(x m , dx ′ m ′ )
where Q is the proposal kernel and α(., .) is the acceptance probability. Two types of transitions are considered here. Adding a marked point to the current configuration is proposed with probability p b and deleting a point from the current configuration is proposed with probability p d . When relevant, the marked object y my to be added to the current marked point configuration x m is chosen with probability density q b (x, y)p M (m y ) with respect to λ Leb X ⊗ λ Leb M . This means that the location of the new object and its mark are simulated independently. Note that the mark is simulated according to the same distribution as that of the marks of the reference process. When relevant, the object in location y to be deleted from current configuration location x m is chosen with probability q d (x, y). We write for (x,

B) in Ξ × F Ξ Q(x m , B) = p b Q b (x m , B) + p d Q d (x m , B).
The so-called "birth" kernel is

Q b (x m , B) = X ×M I(x m ∪ {y my } ∈ B)q b (x, y)p M (m y )dλ Leb (y)dλ Leb (m y ) (10)
Obviously, the new object y my needs to be chosen apart from x m which happens here almost-surely. The so-called "death" kernel is

Q d (x m , B) = y my ∈x m I(x m \ {y my } ∈ B)q d (x, y) (11)
Each proposal move is either accepted of rejected according to an acceptance step detailed later.

If the proposal is accepted, the chain goes to the new state given by the transition kernel. If the proposal is rejected, the chain remains in its current state. Such a transition kernel may be written for A in F Ξ as

P (x m , A) = p b X ×M q b (x, y)p M (m y )α(x m , x m ∪ {y my })I(x m ∪ {y my } ∈ A)dλ Leb X (y)dλ Leb M (m y ) + p d y my ∈x m q d (x, y)α(x m , x m \ {y my })I(x m \ {y my } ∈ A) + I(x m ∈ A) 1 -p b X ×M q b (x, y)p M (m y )α(x m , x m ∪ {y my })dλ Leb X (y)dλ Leb M (m y ) -p d y my ∈x m q d (x, y)α(x m , x m \ {y my })
If the current configuration is the empty configuration, obviously the transition kernel cannot propose to delete a marked point. In this case, the above transition kernel becomes

P (∅, A) = p b X ×M q b (∅, y)p M (m y )α(∅, {y my })I({y my } ∈ A)dλ Leb X (y)dλ Leb M (m y ) + I(∅ ∈ A) 1 -p b X ×M q b (∅, y)p M (m y )α(∅, {y my })dλ Leb X (y)dλ Leb M
In the appendix, we determine the Green ratios associated with our choice of inhomogeneous marked Poisson process as reference. These Green ratios are given by r b (x m , y my ) = h(x m ∪ {y my })p d q d (x ∪ {y}, y)λ(y) h(x m )p b q b (x, y) and

r d (z, y) = h(z \ {y})p b q b (z \ {y}, y) h(x)p d q d (z, y)λ(y) .
Note both the scaling of q b by λ(.) and the absence of p M in these expressions. This comes from the fact that the proposition kernel uses the mark distribution P M of the reference process as mark proposition distribution. An acceptance step then makes use of the Green ratios to determine in a stochastic manner whether or not the proposal should be accepted or rejected in which case the chain remains in its current state. The main steps of the Metropolis-Hastings-Green algorithm are summarized in Algorithm 2.

The definition of p M has already been given in subsection 2.1. We now detail our choice of q d and q b . The death is proposed uniformly among the current configuration, which is a classical choice. This yields

q d (x, y) = 1 n(x)
.

The novelty of our algorithm lies in the birth proposal choice. In previous work (see eg Descombes et al. [START_REF] Descombes | Stochastic geometry for image analysis[END_REF]), birth is proposed uniformly in X . But it would make more sense to favour birth where an Algorithm 2 Metropolis Hastings Green Propose y ∼ q b (x, .)

5:

Propose m y ∼ p M (.)

6:

Compute Green ratio r b (x m , y my )

7:

Accept y my with probability α(x m , x m ∪ {y my }) = min(1, r b (x m , y my )). 8: else 9:

Propose y ∼ q d (x, .) object is more likely to be. Since it is possible to get approximate information about the location of the objects in the image by image processing techniques, it seems logical to use this information for that purpose. In the cell case, we propose to use the distribution given by the normalized circular Hough transform obtained after taking the local maxima over a prespecified range of radii. Doing this may yield a distribution putting zero weight over quite large areas of the image domain, which may be seen quite often in practice since cell images may be very poorly contrasted. In this case, we advocate for rather using a finite mixture distribution having as two components the preceding normalized circular Hough transform (with an important weight eg 0.8) and the uniform distribution in X (with a minor weight eg 0.2) so that the normalized circular Hough transform dominates but the whole space is explored anyway. An example of point sampling according to this proposal is displayed in figure 3. In the fiber case, we use the normalized Radon transform. An example of point sampling according to this proposal is displayed in figure 4.

Here appears another benefit of our choice of inhomogeneous marked Poisson process as reference. Descombes et al. [START_REF] Descombes | Stochastic geometry for image analysis[END_REF] mentioned the difficulty that proposing attractive configurations was not enough to improve the convergence of the algorithm. This was due to the fact that attractive configuration have small Green ratio and the proposition was therefore very likely to be rejected. We recall that in Descombes et al. [START_REF] Descombes | Stochastic geometry for image analysis[END_REF], the reference Poisson process is taken as homogeneous with constant intensity equal to λ > 0. In this case, when death proposal is uniform in X and when birth proposal is generated according to a distribution having density denoted as p(.) with respect to λ Leb X , the birth Green ratio writes down as

r b (x m , y) = h(x m ∪ {y my }) h(x m ) λ (n(x) + 1)p(y)
.

The larger p(y) is, the more likely there is an object in the image at location y. But the larger p(y) is, the smaller r b is. Though counter-intuitive at a first look, this entails that any reasonable proposition is likely to be rejected. To emphasize the difference with our algorithm, we recall that our reference Poisson process is inhomogeneous and consider the case where its intensity function is proportional to the birth distribution. We then get

r b (x, y) ∝ h(x m ∪ {y my }) h(x m ) 1 (n(x) + 1)
so that the intensity function compensates for the birth proposal density term and prevents the birth Green ratio from becoming too small when the proposal is very likely according to p(.).

Results and discussion

In this section, we present the results obtained with the proposed methodology. In the cell case, we compare these results with the results of the reference algorithm that uses uniform birth and death kernels together with a unit rate homogeneous Poisson process as reference process. In the fiber case, we were unable to provide a similar comparison due to the lack of existing algorithm for thick straight line objects. We analyse the results in terms of visual quality of object detection and in terms of computation time and convergence to stability of our algorithm. We also provide the required morphology characterisation. The results presented hereafter are obtained with a dual-core laptop having 2.10 GHz CPU twice.

The parameters are calibrated by basic trial and error method since existing estimation methods did not give satisfying results. The issue of estimating the model parameters is left for future research. For simplicity, we take throughout p b = p d = 1 2 . In the circle case, the parameters of the area interaction process are β = 1, γ d = 30, γ p = 50. The initial configuration is the empty configuration. The parameters specific to the simulated annealing part are T 0 = 0.1 as initial temperature, c = 0.9999 as multiplicative factor in the geometric cooling schedule and N = 50000 as prespecified iteration number. The circular Hough transform needed for our sampling step was computed within the radius range [r min , r max ] with r min = 10 and r max = 25 and then normalised to get a proper distribution. The location sampling distribution was then taken as the finite mixture having as two components the normalised circular Hough transform and the uniform distribution with respective weight 0.8 and 0.2. We ran the reference algorithm with the uniform distribution on [r min , r max ] as radius distribution with r min = 10 and r max = 25. The data term parameters were taken as ρ = 15 and d 0 = 40.

When implementing our method in the fiber case, the parameters of the area interaction process are taken as β = 1, γ d = 50 and γ p = 1000. The initial configuration is the empty configuration. The parameters specific to the simulated annealing part are T 0 = 10 as initial temperature, c = 0.9999 as multiplicative factor in the cooling schedule and N = 100000 as prespecified iteration number. The fiber width was simulated uniformly on [w min , w max ] with w min = 15 and w max = 30. The data term parameters were taken as s = 1, s ′ = 5000, V max = 1 and τ = 10.

The detected objects are represented on the original image on figure 5 in the cell case and on figure 8 in the fiber case. A visual examination leads to the conclusion that our method provides a valuable object detection.

For the cell image, our algorithm runs in 43 minutes for 50000 iterations while the reference algorithm runs in 41 minutes with the same iteration number making the computation times similar. For the fiber image, our algorithm runs in 18 minutes for 100000 iterations.

The evolution of the number of detected cells is represented as a function of the number of iterations in the cell case on figure 6 and in the fiber case on figure 9. Noticeably, in the cell case, our algorithm converges quickly to equilibrium in a stable way which eases the convergence diagnosis.

The required morphology characterisation is provided in figure 7 in the cell case and in figure in the fiber case 10. Noticeably, in the cell case, our results illustrate the fact that the estimated radius distribution shares a similar shape with the mark proposal distribution. This highlights the importance of calibrating the proposal in a data-driven way when possible.

As a conclusion, our results illustrate the reliability and the efficiency of our data-driven technique to guide the Markov chain search and get quickly the required morphology characterisation.

Appendix: some mathematical details

We present here some technical results to prove the convergence of the proposed algorithm. We follow the lines of Stoica [START_REF] Stoica | Marked point processes for statistical and morphological analysis of astronomical data[END_REF].

More generally, we show the convergence of the proposed algorithm when the following conditions are simultaneously fulfilled:

• the target distribution π is a locally stable process having density h(.) with respect to an inhomogeneous Poisson process with intensity function λ(.),

• the birth and death probabilities satisfy 0 < p b , p d < 1 with 0 < p b + p d < 1,

•

u n := sup y my ∈X ×M,x m ∈Ξn q d (x, y)λ(y) q b (x \ {y}, y) n→∞ -→ 0 , (12)
and q b (x \ {y}, y) q d (x, y)

1 λ(y) > 0, ∀x m ∈ Ξ, ∀y my ∈ X × M . (13)
First, the algorithm should reproduce a transition kernel P such that the target distribution π is an invariant distribution of the chain. It turns out that the target distribution π is an invariant distribution of the chain if the proposal kernel P and the acceptance probabilities α(., .) are chosen so that the Markov chain is reversible with respect to π. The reversibility property is fulfilled when for any

A, B ∈ F Ξ πP (A, B) = A B I A (x m )I B (x ′ m ′ )π(dx m )P (x m , dx ′ m ′ )
is symmetric when exchanging the sets A and B. We require that the proposal kernel and the acceptance probabilities satisfy the stronger property of detailed balance condition which writes down here as

p b h(x m )q b (x, y)α(x m , x m ∪ {y my }) = p d h(x m ∪ {y my })q d (x ∪ {y}, y)α(x m ∪ {y my }, x m )λ(y) (14)
where λ(.) is the intensity function of the reference inhomogeneous Poisson process. We first show that the detailed balance condition [START_REF] Lavielle | Controlled formation of polycaprolactone ultrathin electrospun nanofibers by hydrolytic degradation approach[END_REF] implies the reversibility of the transition kernel.

Since the proposal kernel Q only enables moves from Ξ n to Ξ n+1 and from Ξ n+1 to Ξ n and since the acceptance step only allows in addition the possibility of no move, we only need to prove that πP (A n , B n+1 ) = πP (B n+1 , A n ) for A n ⊆ Ξ n and B n+1 ⊆ Ξ n+1 for any n. On the one hand, the definition of πP gives

πP (A n , B n+1 ) = (X ×M) n I(x m ∈ A n )P (x m , B n+1 )π(dx m ) = x m ∈An x ′m ′ ∈B n+1 α(x m , x ′ m ′ )p b Q b (x, dx ′ m ′ )π(dx m )
where π(dx m ) = h(x m )π Λ (dx m ) and where Q b is given by [START_REF] Konomi | Bayesian object classification of gold nanoparticles[END_REF]. We get

πP (A n , B n+1 ) = (X ×M) n+1 I(x m ∈ A n , x m ∪ {y my } ∈ B n+1 ) × α(x m , x m ∪ {y my })p b q b (x, y)p M (m y )h(x m )λ Leb X (dy)λ Leb M (dm y )π Λ (dx m )
Now, we substitute π Λ by its expression in (2). This gives

πP (A n , B n+1 ) = exp(-Λ(X )) n! (X ×M ) n+1 I(x m ∈ A n , x m ∪ {y my } ∈ B n+1 )α(x m , x m ∪ {y my }) × p b q b (x, y)p M (m y )h(x m )λ Leb X (dy)λ Leb M (dm y ) n i=1 λ(x i )p M (m i )λ Leb X (dx i )λ Leb M (dm i )
On the other hand, we have

πP (B n+1 , A n ) = (X ×M) n+1 I(z m ∈ B n+1 )P (z m ∈ A n )π(dz m ) = z m ∈B n+1 z ′m ′ ∈An α(z m , z ′ m ′ )p d Q d (z m , dz ′ m ′ )π(dz m )
where π(dz m ) = h(z m )π Λ (dz m ) and where Q d is given by [START_REF] Lacoste | A Comparative Study of Point Processes for Line Network Extraction in Remote Sensing[END_REF]. We get

πP (B n+1 , A n ) = (X ×M) n+1 y my ∈X ×M I(y my ∈ z m )I(z m ∈ B n+1 , z m \ {y my } ∈ A n ) × α(z m , z m \ {y my })p d q d (z, y)h(z m )π Λ (dz m )
We reverse the order of the sum and the integral, all the terms being positive:

πP (B n+1 , A n ) = y my ∈X ×M (X ×M) n+1 I(y my ∈ z m )I(z m ∈ B n+1 , z m \ {y my } ∈ A n ) × α(z m , z m \ {y my })p d q d (z, y)h(z m )π Λ (dz m )
Now, we substitute π Λ by its expression on (2). This gives

πP (B n+1 , A n ) = exp(-Λ(X )) (n + 1)! y my ∈X ×M (X ×M) n+1 I(y my ∈ z m )I(z m ∈ B n+1 , z m \ {y my } ∈ A n ) × α(z m , z m \ {y my })p d q d (z, y)h(z m ) n+1 i=1 λ(z i )dλ Leb X (z i )p M (dm i )dλ Leb M (m i ) Letting x m = z m \ {y my } in the integrals, we get πP (B n+1 , A n ) = exp(-Λ(X )) (n + 1)! y my ∈X ×M (X ×M) n+1 I(x m ∪ {y my } ∈ B n+1 , x m ∈ A n )α(x m ∪ {y my }, x m ) × p d q d (x ∪ {y}, y)h(x m ∪ {y my }) n i=1 λ(x i )dλ Leb X (x i )p M (dm i )dλ Leb M (m i )λ(y)dλ Leb X (y)p M (dm y )dλ Leb M (m y )
Only the terms for which the indicator I(x m ∪ {y my } ∈ B n+1 , x m ∈ A n ) is non zero are non-zeros. There are (n + 1) such terms in which case they are all equal to the same integral (since it does not depend on the integration order). This entails

πP (B n+1 , A n ) = exp(-Λ(X )) n! (X ×M) n+1 I(x m ∪ {y my } ∈ B n+1 , x m ∈ A n )α(x m ∪ {y my }, x m ) × p d q d (x ∪ {y}, y)h(x m ∪ {y my }) n i=1 λ(x i )dλ Leb X (x i )p M (dm i )dλ Leb M (m i )λ(y)dλ Leb X (y)p M (dm y )dλ Leb M (m y )
The detailed balance condition in [START_REF] Lavielle | Controlled formation of polycaprolactone ultrathin electrospun nanofibers by hydrolytic degradation approach[END_REF] gives the reversility which concludes this proof.

We determine the Green ratios for both the birth and death proposals according to the method of Green [START_REF] Green | Reversible jump Markov Chain Monte Carlo computation and Bayesian model determination[END_REF]. To this aim, we define a symmetric measure ψ : F Ξ × F Ξ → R + such that the Radon-Nikodym derivative of πQ with respect to ψ exists. This derivative is denoted by D. The proof runs along similar lines than the above one but the calculus are presented for the sake of clarity and completeness. Let us recall that the measure πQ is defined for (A, B) in

F Ξ × F Ξ by πQ(A, B) = x m ∈A Q(x m , B)π(dx m ).
Since the proposal kernel Q only enables moves from Ξ n to Ξ n+1 and from Ξ n+1 to Ξ n , we only have to study the derivability of πQ on sets (A n , B n+1 ) or (A n+1 , B n ) with A n , B n ⊆ Ξ n and A n+1 , B n+1 ⊆ Ξ n+1 for any arbitrary n. On the first hand,

πQ(A n , B n+1 ) = (X ×M) n I(x m ∈ A n )Q(x m , B n+1 )π(dx m ) = (X ×M) n I(x m ∈ A n )p b Q b (x m , B n+1 )π(dx m )
with π(dx m ) = h(x m )π Λ (dx m ) and where Q b is given by [START_REF] Konomi | Bayesian object classification of gold nanoparticles[END_REF]. We get

πQ(A n , B n+1 ) = (X ×M) n+1 I(x m ∈ A n , x m ∪ {y my } ∈ B n+1 ) p b q b (x, y)h(x m )λ Leb X (dy)p M (m y )λ Leb M (dm y )π Λ (dx m )
We define a measure ψ + n on Ξ n × Ξ n+1 by

ψ + n (A n , B n+1 ) = (X ×M) n+1 I(x m ∈ A n , x m ∪ {y my } ∈ B n+1 ) λ(y)λ Leb X (dy)p M (m y )λ Leb M (dm y )π Λ (dx m ) D b (x m , x m ∪ {y my }) = p b q b (x, y)h(x m ) λ(y)
On the other hand,

πQ(A n+1 , B n ) = (X ×M) n+1 I(z m ∈ A n+1 )Q(z m , B n )π(dz m ) = (X ×M) n+1 I(z m ∈ A n+1 )p d Q d (z m , B n )π(dz m )
with π(dz m ) = h(z m )π Λ (dz m ) and where Q d is given by [START_REF] Lacoste | A Comparative Study of Point Processes for Line Network Extraction in Remote Sensing[END_REF]. We get

πQ(A n+1 , B n ) = (X ×M) n+1 y my ∈z m I(z m ∈ A n+1 , z m \ {y my } ∈ B n )p d q d (z, y)h(z m )π Λ (dz m )
We reverse the order of the sum and the integral, all the terms being positive:

πQ(A n+1 , B n ) = y my ∈X ×M (X ×M) n+1 I(z m ∈ A n+1 , y my ∈ z m , z m \{y my } ∈ B n )p d q d (z, y)h(z m )π Λ (dz m )
Letting x m = z m \ {y my } in the integrals, we get

πQ(A n+1 , B n ) = y my ∈X ×M (X ×M) n+1 I(x m ∪ {y my } ∈ A n+1 , x m ∈ B n ) p d q d (x ∪ {y}, y)h(x m ∪ {y my })λ(y)λ Leb X (dy)p M (m y )λ Leb M (dm y )π Λ (dx m )
since in view of (2), we have that π Λ (dz m ) = π Λ (dx m )λ(y)λ Leb X (dy)p M (m y )λ Leb M (dm y ). Only the terms for which the indicator I(x m ∪ {y my } ∈ A n+1 , x m ∈ B n ) is non zero are non-zeros. There are (n + 1) such terms in which case they are all equal to the same integral (since it does not depend on the integration order). This entails

πQ(A n+1 , B n ) = (n + 1) (X ×M) n+1 I(x m ∪ {y my } ∈ A n+1 , x m ∈ B n ) p d q d (x ∪ {y}, y)h(x m ∪ {y my })λ(y)λ Leb X (dy)p M (m y )λ Leb M (dm y )π Λ (dx m )
We define a measure ψ - n on Ξ n+1 × Ξ n by

ψ - n (A n+1 , B n ) = (n + 1) (X ×M) n+1 I(x m ∪ {y my } ∈ A n+1 , x m ∈ B n ) λ(y)λ Leb X (dy)p M (m y )λ Leb M (dm y )π Λ (dx m )
The Radon-Nikodym derivative of πQ with respect to ψ - n on Ξ n+1 × Ξ n is then defined by

D d (x m ∪ {y my }, x m ) = p d q d (z, y)h(x m ∪ {y my })
We then define a measure ψ on Ξ × Ξ as the measure concentrated on the set

∞ n=0 disjoint ((Ξ n × Ξ n+1 ) ∪ disjoint (Ξ n+1 × Ξ n ))
which coincides with ψ + n on Ξ n × Ξ n+1 and with ψ - n on Ξ n+1 × Ξ n . It is straightforward to see that this actually defines a measure on Ξ × Ξ. To prove the symmetry of ψ, we only have to prove that

ψ(A n+1 , B n ) = ψ(B n , A n+1 ) for A n ∈ Ξ n and B n+1 ∈ Ξ n+1 . This is immediate since ψ(B n , A n+1 ) = ψ + n (B n , A n+1 ) = exp(-Λ(X )) n! (X ×M) n+1 I(x m ∈ B n , , x m ∪ {y my } ∈ A n+1 )λ(y)λ Leb X (dy) × p M (m y )λ Leb M (dm y ) n i=1 λ(x i )λ Leb X (dx i )p M (m i )λ Leb M (dm i ) and ψ(A n+1 , B n ) = ψ - n (A n+1 , B n ) = (n + 1) exp(-Λ(X )) (n + 1)! (X ×M) n+1 I(x m ∪ {y my } ∈ A n+1 , x m ∈ B n ) × λ(y)λ Leb X (dy)p M (m y )λ Leb M (dm y ) n i=1 λ(x i )λ Leb X (dx i )p M (m i )λ Leb M (dm i )
The derivative D exists and coincides with D b on Ξ n × Ξ n+1 and with D d on Ξ n+1 × Ξ n . We then deduce the following Green ratios:

r b (x m , y) = D d (x m ∪ {y my }, x m ) D b (x m , x m ∪ {y my }) = h(x m ∪ {y my }) h(x m ) p d p b q d (x ∪ {y}, y) q b (x, y) λ(y) and r d (z m , y my ) = D d (z m \ {y my }, z m ) D b (z m , z m \ {y my }) = h(z m \ {y my }) h(z m ) p b p d q b (z \ {y}, y) q d (z, y) 1 λ(y)
The acceptance probabilities are then defined in the birth case as

It is easy to check that r b (x m , x m ∪ {y my }) = 1 r d (x m ∪ {y my }, x m ) so that the detailed balance condition in ( 14) is satisfied.

A potential problem is that the condition of reversibility alone implies the initial distribution needs to be chosen according to the distribution of interest. This problem is solved if the chain is irreducible, aperiodic and Harris recurrent.

Irreducibility is a crucial property in the set-up of MCMC algorithms. Roughly speaking, a φirreducible Markov chain is able to reach any set A which is large enough in the sense that φ(A) > 0. For any positive integer n, let P n be the n-step transition probability kernel which is defined inductively for x m ∈ Ξ and B ∈ F Ξ by

P n (x m , B) = Ξ P (x m , dz m )P n-1 (z m , B)
where P 0 (x m , B) = δ B (x m ) with δ B being the Dirac measure concentrated at B. We show that our chain satisfies the irreducibility property provided the conditions ( 13) and ( 12) are fulfilled. To this aim, we exhibit a positive measure φ such that, for any set A in F Ξ such that φ(A) > 0, there exists a positive integer n 0 such that P n 0 (x m , A) > 0 for any x m in Ξ. We choose as positive measure the Dirac measure concentrated at ∅, that is φ = δ {∅} , so that the only set A in F Ξ having positive measure is {∅}. The Markov property of the chain entails that, for any x m ∈ Ξ, for any positive integer n 0 such that n 0 ≥ n(x), we have P n 0 (x m , {∅}) ≥ P n(x) (x m , {∅})P n 0 -n(x) (∅, {∅}) where

P (∅, {∅}) = 1 -p b X ×M q p (∅, y)p M (m y )α(∅, {y my })dλ Leb X (y)dλ Leb M (m y ) ≥ 1 -p b = p d
and where

P n(x) (x, {∅}) ≥ p n(x) d min z m ∈Ξ,y my ∈z m : n(z)≤n 0 {α(z m , z m \ {y my })} n(x)
.

For z m in Ξ with n(z) ≤ n 0 and y my ∈ z m , the death acceptance probability is defined by [START_REF] Meyn | Markov chains and stochastic stability[END_REF]. Since the target distribution is a locally stable process, there exists a K > 0 such that for z m in Ξ and

y my ∈ z m , h(z m \ {y my }) h(z m ) ≥ 1 K .
On the other hand, letting

∆(n 0 ) = min z m ∈Ξ,y my ∈z m : n(z)≤n 0 q b (z \ {y}, y) q d (z, y) 1 λ(y) ,
we deduce that for any x m in Ξ such that n(x) ≤ n 0 , for any

y my in x m α(x m , x m \ {y my }) ≥ min 1, ∆(n 0 ) K p b p d
In view of [START_REF] Lacoste | Road network extraction in remote sensing by a Markov Object process[END_REF], one may find a n 0 large enough so that

∆(n 0 ) < 1 K p b p d .
in which case for any x m in Ξ such that n(x) ≤ n 0 , for any

y my in x m α(x m , x m \ {y my }) ≥ ∆(n 0 ) K p b p d .
We get that for any x m ∈ Ξ, one may find a n 0 large enough so that

P n 0 (x m , {∅}) ≥ p n(x) d ∆(n 0 ) K p b p d n 0 p n 0 -n(x) d = p b ∆(n 0 ) K n 0 . (16)
In view of [START_REF] Lavielle | Structuring and molding of electrospun nanofibers: Effect of electrical and topographical local properties of micro-patterned collectors[END_REF], this gives the φ-irreducibility of P .

The Harris recurrence is the property that guarantees the non-existence of sets of configurations Ξ ′ ⊂ Ξ such that π(Ξ ′ ) = 1 and lim n→∞ sup A∈F Ξ |P n (x m , A)π(A)| = 0 for each x m in Ξ ′ . Hence the recurrence property guarantees the asymptotic convergence of a MCMC sampling algorithm regardless of the initial conditions. According to the paper of Meyn and Tweedie [START_REF] Meyn | Markov chains and stochastic stability[END_REF], a sufficient condition for Harris recurrence is the drift condition. The drift condition says that, if the transition kernel P is φ-irreducible, if there is a small set C in F Ξ and a function unbounded off small sets V : X i → R + such that Ξ P (x m , dz m )V (z m ) ≤ V (x m ), ∀x m / ∈ C then P is Harris recurrent. We recall that a set C in F Ξ is small if there are a positive integer n 0 and a non-zero measure ν such that

P n 0 (x m , B) ≥ ν(B), ∀x m ∈ C, ∀B ∈ F Ξ
First, set a number A > 1 and then set a number 0 < ε < min 1, p d p b A . In view of [START_REF] Lacoste | Road network extraction in remote sensing by a Markov Object process[END_REF], one may find a N ε large enough so that for ny n ≥ N ε

p d p b Ku n ≤ ε (17)
Equation [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF] shows that the set

C(ε) = {x m ∈ Ξ : n(x) ≤ N ε }
is small with the choice ν(.) = (p d ∆(N ε )/K) Nε δ {∅} (.). We also recall that a function V : Ξ → R + is unbounded off small sets if the level sets

≤ V (x m ) .

So far, we know that the Harris recurrent chain should reach the equilibrium asymptotically but we do not know when this will happen. One way of solving this problem is to build an ergodic chain (that is Harris recurrent and aperiodic). In this case, there is no need to wait infinitely until the equilibrium. The property of aperiodicity is immediate since the self-transition probabilities The last refinement consists of proving that the geometric drift property is fulfilled so that the convergence of the chain to its equilibrium distribution happens at a geometric speed. The geometric drift condition says that a φ-irreducible and aperiodic Markov chain with transition kernel P is geometrically ergodic if there exists a function V : Ξ → (1, ∞), two numbers b in R and a < 1 and a small set C in F Ξ such that Ξ V (z m )P (x m , dz m ) ≤ aV (x m ) + bI C (x m ), ∀x m ∈ Ξ .

P (x m , {x m }) = 1 -p b X ×M q p (x,
Taking the set C(ε) and the function V as previously defined for some 1 < A < 1 + 1/p b and taking b = A Nε+1 gives the desired property in view of [START_REF] Robert | Monte Carlo Statistical Methods[END_REF].

Indeed, if x m / ∈ C(ε), then (18) reads as ≤ p b (A -1)

:=a ′′ <1

A n(x) + A Nε+1

Taking a = max(a ′ , a ′′ ) gives the result.

Figure 1 :

 1 Figure 1: A Scanning Electron Microscopy image of polymer fibers and an optical microscopy image of cells. Image courtesy by Guy Schlatter, ICPEES, University of Strasbourg, France.

.

Algorithm 1 annealing 1 :

 11 Simulated Choose an initial temperature T = T 0 2: while maximal number of iterations not reached do 3:

10 : 11 :

 1011 Compute Green ratio r d (x m , y my ) Accept y my with probability α(x m , x m \ {y my }) = min(1, r d (x m , y my )).

Figure 3 :

 3 Figure 3: A non-uniform sampling example in the cell case.

  Its Radon transform (c) A sampling example of 100000 points according to the normalized Radon transform

Figure 4 :

 4 Figure 4: A non-uniform sampling example in the fiber case.

  α(x m , x m ∪ {y my }) = min(1, r b (x m , y my )) and in the death case as α(z m , z m \ {y my }) = min(1, r d (z m , y my )) .

C 1 y 1 . 1 A - 1 + 1 V

 11111 V (α) = {x m ∈ Ξ : V (x m ) ≤ α} are small for any α > 0. The choice V (x m ) = A n(x) entails as previously that the setsC V (α) = {x m ∈ Ξ : n(x) ≤ log α A }are small for any α > 0. It remains to get the negativity of the drift. Applying the transition kernel P to the function V leads toΞ V (z m )P (x m , dz m ) = p b A n(x) (A -1) X ×M q b (x, y)α(x m , x m ∪ {y my })p M (m y )dλ Leb X (y)dλ Leb M (m y ) + p d A n(x) 1 A my ∈x m q d (x, y)α(x m , x m \ {y my }) + A n(x)In view of (17), we deduce for x m / ∈ C(ε) that α(x m , x m ∪ {y my }) ≤ ε and that α(x m , x m \ {y my }) = This gives for x m / ∈ C(ε) Ξ V (z m )P (x m , dz m ) ≤ p b (A -1)ε + p d (x m )

  y)p M (m y )α(x m , x m ∪ {y my })dλ Leb X (y)dλ Leb M (m y ) p d y my ∈x m q d (x, y)α(x m , x m \ {y my })are positive for x m in A such that π(A) > 0.

ΞV 1 A - 1 + 1 : 1 y

 1111 (z m )P (x m , dz m ) ≤ p b (A -1)ε + p d =a ′ <1 V (x m ) + b × 0 . If x m ∈ C(ε), then Ξ V (z m )P (x m , dz m ) = p b A n(x) (A -1) X ×M q b (x, y)α(x m , x m ∪ {y my })p M (m y )dλ Leb X (y)dλ Leb M (m y ) ≤1 + p d A n(x) 1 A my ∈x m q d (x, y)α(x m , x m \ {y my }) <0 +A n(x)

( a )

 a Detected cells using our method (b) Detected cells using the reference algorithm

Figure 5 :Figure 6 :

 56 Figure 5: Detected cells in the same SEM image

( a )

 a Distribution used to simulate cell radius using our method (b) Estimated distribution of cell radius using our method (c) Uniform distribution used to simulate cell radius using the reference algorithm (d) Estimated distribution of cell radius using the reference algorithm

Figure 7 :Figure 8 :

 78 Figure 7: Histograms of cell radius distribution

Figure 9 :

 9 Figure 9: Number of detected fibers as a function of iteration number using our method

Figure 10 :

 10 Figure 10: Histogram of fiber width and polar histogram of fiber orientation using our method

1 :

 1 Let x m be the current configuration 2: Propose birth with probability p b and death with probability p d 3: if birth then

	4:

Acknowledgements: The first author is on a post-doctoral position financed by the Agence Nationale de la Recherche (ANR) via the project ANR 2011 NANO 018 01 NeoTissage. The second author also benefits from financial support from the project ANR 2011 NANO 018 01 NeoTissage.


Martien Den Hertog
email: martien.den-hertog@grenoble.cnrs.fr

Rudeesun Songmuang

Eva Monroy

Polarization fields in GaN/AlN nanowire heterostructures studied by Off axis holography

In this work, we present an off-axis holography study of GaN/AlN heterostructured nanowires grown by plasma-assisted molecular-beam epitaxy. We discuss the sample preparation of nanowire samples for electron holography and combine potential profiles obtained using holography with theoretical calculations of the projected potential in order to gain understanding of the potential distribution in these nanostructures. The effects of surface states are discussed.

Introduction

Semiconductor nanowires (NWs) have become a powerful and diverse family of functional materials for nano-electronics, optoelectronics and biotechnology. III-nitride semiconductors, with their direct band gap, flexibility for band engineering, and capability of operating in extreme environments, are promising materials for single-NW devices, namely high-resolution chemical sensing or high-spectralcontrast photodetectors. A particularity of nitride materials is the very strong polarization-related internal electric field present in heterostructures, which can be assessed through its effect on optical or electrical properties. A direct measurement of the internal electric field profile is more challenging but could in principle be realized by off-axis electron holography. Holography on nitride materials was already reported in [1,2,3,4]. In addition to the technical difficulties to measure the field distribution, these direct measurements are necessarily perturbed by the superposition of information on surface depletion due to surface states and by the effect of strain on the internal field distribution. Taking these phenomena into account is critical to understand the physical properties of such nanostructures.

In this paper we use medium resolution off-axis electron holography to measure the projected potential in GaN/AlN NW heterostructures and compare our experimental results with threedimensional (3D) potential simulations to discriminate the 3D field distribution and the effect of surface charge. Furthermore, we describe the challenges of holography on (nitride) NW heterostructures, which have to be faced to allow quantitative mapping of the potential.

Off-axis medium resolution electron holography was performed on an FEI Titan 80-300 working at 200 kV. Plasma-assisted molecular beam epitaxy (PAMBE) is a well-established technique to synthesize nanometer-size catalyst-free GaN NWs growing along the c crystallographic axis [5,6]. This growth method allows n-type and p-type doping [7], as well as band engineering by heterostructuring [8,9,10,11,12]. For this study GaN NWs containing a 10-period GaN/AlN (7 nm / 7 nm) superlattice, followed by 30 nm of GaN and 30 nm of AlN, were grown on Si(111) patterned pillars by PAMBE (see Fig. 1). Growth proceeded under N-rich conditions (Ga/N ratio ~ 0.2-0.4) at a substrate temperature ~ 790 ºC. An AlN buffer layer is used to improve the vertical alignment of the NWs [13]. The NWs under consideration have a length of 1.2-1.5 μm and a diameter of 30-80 nm. The sample for electron holography was prepared using the cleavage technique [14,15], which allowed studying epitaxial NWs while they were still attached to their substrate and suspended in vacuum, without any modification (electrical or structural) of the NW.

The sample preparation is crucial to be able to obtain good holography data from NW samples. In samples with a high density of NWs, the post geometry has several advantages over conventional dispersed NW samples on a carbon grid. First of all, a very precise orientation of the NW several degrees off the zone axis, to minimize diffraction contrast, is possible as the NWs are grown epitaxially on the Si post. Second, it is very easy to find back the same single NW as the posts also serve as markers. It should be noted that for lower density NW systems, a cleaved sample of NWs on a planar substrate can be sufficient to study free-standing NWs without overlap with other NWs close to the cleaved edge [14,15]. However, as GaN NWs grown by PAMBE are very closely packed, we need to find NWs growing at the edge of the post to prevent superposition of other NWs. Previously we have demonstrated that GaN NWs are hexahedral structures with relatively flat {10-10} m-plane sidewall facets, and they exhibit N-polarity, as identified by annular bright field (ABF-) and high angle annular dark field scanning transmission electron microscopy (HAADF-STEM) [16]. As off-axis holography is sensitive to the projected potential we expect to observe the effect of the polarity on the potential in the structure. Figure 2(a) presents a high angle annular dark field (HAADF) STEM image of a GaN/AlN heterostructured NW. The axial superlattice is enveloped by an AlN shell arising from the lateral growth during the deposition of the AlN sections. As a result of the elevated growth temperature, a certain interdiffusion is observed at the interfaces. This phenomenon is particularly important when starting the growth of AlN on GaN, giving rise to a graded interface, which extends along about 1.5 nm. Figure 2b shows the phase image obtained by off-axis electron holography, tilted several degrees off the Si [112] zone axis (ZA). The fringe spacing was 1 nm, giving a spatial resolution in the reconstructed phase image of around 3 nm. The various GaN and AlN regions can be clearly resolved. In the absence of magnetic fields and diffraction contrast, the phase difference Δφ between electrons that have passed through vacuum and electrons that have traversed the specimen is related to the potential V(x, y, z), by

Δϕ =C E V x, y,z ( )dz 0 t ∫ , (1)
where C E is an interaction constant which depends on the acceleration voltage (7.29×10 6 rad V -1 m -1 at 200 kV), and t is the thickness of the sample in the beam direction z. V(x, y, z) is a combination of the mean inner potential (MIP) of the pure components, namely GaN and AlN (MIP = 16.8 V and 14.2 V, respectively [17]), and of the electrostatic potential, related to the electric field, reflecting the combined effects of spontaneous and piezoelectric polarization and surface charges. In Fig 2(c), a projected potential profile is obtained from the phase image in Fig. 2(b) along the white vertical arrow using eq. 1 and the measured NW thickness.

The experimental potential profile is compared to two simulated potential profiles obtained from projected 3D simulations of a NW structure with the layer sequence of the sample in Fig. 2(a) using the nextnano3 software [18] with the material parameters in ref. [19], one including a negative surface charge density of -2x10 13 cm -2 , which is one order of magnitude higher than a value from literature [20], and one without surface charge. These simulations take into account the NW hexagonal cross section with m sidewall facets, with the dimensions extracted from our STEM measurements, including the presence of the AlN shell. The strain distribution is calculated to minimize the elastic energy in the system considering a free-standing NW in vacuum. The calculation of the 3D potential profile takes into consideration the local variations of the band gap due to strain, as well as the 3D internal electric field induced by spontaneous and piezoelectric polarization. However, our simulations do not include the variation in MIP between GaN and AlN; thus, the experimentally obtained potential profile was shifted along the y-axis to allow comparison to the simulated profiles. We verified that the step in phase over the AlN top region is consistent with the thickness and the MIP of AlN. For the studied NW, the expected change in potential when going from a GaN to AlN region due only to the variation in MIP would be around 0.8 V, if we assume the NW diameter of 44 nm is roughly equal to the thickness. The experimentally measured potential change at the GaN/AlN heterointerfaces is much larger than the calculated 0.8 V, indicating that indeed the measurement is sensitive to the electric field. However, in the studied NWs there are composition gradients at the interfaces, creating also a gradient in the MIP. Therefore, we cannot provide a quantitative assessment of the interface charge in this structure, but quantitative measurements are feasible in structures with abrupt interfaces and homogeneous alloy composition.

Making a qualitative analysis, the experimental profile is comparable to the simulated profiles and consistent with the N-polarity of the NW. Interestingly, the internal field is so strong that the surface charge does not have much influence on the projected potential, even though the value of the surface charge used in the simulation is higher than experimentally measured [20]. On the other hand, the measurement at the AlN top region presents a deviation from simulations, which might be due to charging of the NW under the electron beam. To investigate charging issues, we display in Fig. 2d a phase profile in the vacuum at the side of the NW along the black horizontal dotted arrow in Fig 2b .  Going away from the NW the phase decreases, indicating a positive charging of the NW [21]. This could explain the lowering of the potential in the AlN top region.

In summary we have shown first results of holography on GaN/AlN NW heterostructures, compared to theoretical calculations. We confirm the N-polarity of the NW and observe a positive charging of the sample under the electron beam. Comparison with simulated profiles does not allow setting a value on the surface charge of the NW. We have described some of the challenges that have to be met to allow quantitative field mapping in nitride NW heterostructures.

Financial support from the French FMN-SMINGUE 2011, the French CNRS and CEA METSA network, the

Figure 1 .

 1 Figure 1. (a,b,c) Scanning electron microscopy (SEM) top view at increasing magnification of GaN/AlN NW heterostructures grown on Si (111) patterned pillars.

Figure 2 .

 2 Figure 2. (a) HAADF STEM image of a GaN/AlN NW heterostructure. (b) Phase image of the same NW. (c) Potential profile obtained from the phase along the vertical white arrow compared with two simulated profiles one including a surface charge density of -2x10 13 electron charges cm -2 and one without surface charges. (d) Phase profile in the vacuum obtained along the horizontal black dotted arrow.

ANR-2011-NANO-027 "UVLamp" project, the ANR 12 JS10 00201 "COSMOS" project, the EU ERC-StG "TERAGAN" (#278428) project, and the Société Française des Microscopies (SFμ) for a travel grant is acknowledged. We thank Adeline Grenier for providing the Si(111) patterned pillars as the growth substrate.