In: Paparoditis N., Pierrot-Deseilligny M.. Mallet C.. Tournaire O. (Eds). IAPRS. Vol. XXXVIII. Part ЗА - Saint-Mandé, France. September 1-3, 2010 
Algorithmically, the computation of the global minimum of this 
energy is performed by a simulated annealing scheme (Azencott, 
1992). 
2.1.1 Data energy Two main approaches for computing the 
data energy term were previously introduced to extract surface 
networks (Chatelain et al., 2009a). The first one was carried out 
in a Bayesian framework. Nevertheless, it doesn’t sufficiently 
take into account the morphology of objects. The second one is 
more efficient. It consists in defining a local energy U#(u) E 
[—1,1] associated with each object u. Then, the data energy as 
sociated with a configuration x is proportional to the sum of all 
the local energies computed on each object it of x: 
Ue (x.y) = 7 d V Ue (it, y) (3) 
u£x 
where the parameter 7> 0 corresponds to the data energy 
weight. Its definition relies on a contrast measure between the 
distribution of the set of pixels belonging to an object u and the 
distribution of pixels belonging to its border J rp {u). The data en 
ergy Uq (it) is built as a qualification of this contrast measure in 
order to promote well-placed objects and penalize bad ones: 
Ue(u) = Q( 
d(u,P p (u)) 
do 
(4) 
where d(u, J- P (u)) is the Bhattacharya distance between the ob 
ject it and its boundary T p (u). Qa 0 : K + ^ [—1,1] is a quality 
function. It attributes a negative value to objects that have a high 
radiometric distance w.r.t. their border (i.e. d{u. T p {u)) is above 
the threshold do) and a positive value otherwise: 
Q(x) = 
1 - æ 1 / 3 
exp(-~y^) - 1 
if x < 1, 
if X > 1. 
(5) 
Using a cubic root in this quality function ensures a moderate pe 
nalization when the distance d(u. J rp {u)) is close to the threshold 
do. To illustrate the behavior of such a function, we give its plot 
in figure 2(left). 
2.1.2 Prior energy This energy corresponds to a penalization 
of the overlapping objects and thus avoids detecting the same ob 
ject several times. For this purpose, we define the neighborhood 
system corresponding to overlapping objects by the following re 
lationship: 
2.2 Estimation algorithm 
Before revealing the proposed estimation method, it is necessary 
to identify the parameters to be estimated. The parameter vec 
tor 9 embedded in the considered model is composed of: the 
activity parameter 0, the maximal overlapping ratio s, the data 
energy weight 7the threshold do involved in the quality func 
tion and the width of the boundary of an object p. Some of these 
parameters can be set in a deterministic way as they have an in 
tuitive interpretation: The boundary width p is usually set to 1 or 
2. the threshold s is set as the tolerated recovery proportion be 
tween objects. The previous study of the model parameters has 
shown that the activity parameter 0 and the weight parameter 7a 
are correlated. The experience proved that a low value of 0 is 
compensated by a high value of 77 and vice versa. Therefore, 
w'e decide to estimate only one parameter which is the weight 7d 
and set manually the other ones. The problem of parameter esti 
mation consists in identifying the most likely parameter vector 9 
that maximizes the image likelihood fe{y). As the density of the 
observation fe{y) is unknown, we propose to maximize the ex 
tended likelihood fe(x, y). However, the configuration x is like 
wise unknown. The study of estimation methods (Chatelain et 
al., 2009a) showed that the stochastic version of the Expectation- 
Maximization (SEM) algorithm (Celeux et al., 1996) is very rele 
vant in such a situation. It consists in iterating the three following 
steps: 
(8) 
(9) 
(10) 
Nevertheless, the E step is unreachable as the extended density 
fo(x ik) . y) is not tractable (the normalizing constant is inacces 
sible). To overcome this difficulty, we propose to approximate the 
extended likelihood by the pseudo-likelihood (Besag, 1975, Bad- 
deley and Turner, 2000). Its expression in the case of incomplete 
data is given by: 
1. S step: Simulate a configuration of objects: 
x [k) ~ f 6 k(X/y) 
2. E Step: Evaluate the extended log-likelihood: 
Q(9,9 k :y) = logfe(x (k \y) 
3. M step: Maximize the log-likelihood: 
9 k + ] = arg max Q{9. 9 k ; y) 
Va’j, Xj E W, x, 
Area(37 П Xj) 
min(Area(37), Area (ж,)) 
< s 
(6) 
The parameter s E [0,1] represents the maximum recovery ratio 
between two objects. Afterwards, we introduce a “Hard Core” 
process w'hich penalizes any configuration having at least one 
overlapping object pair with a recovery rate exceeding a thresh 
old s. The non-normalized density of such a process is written 
as follows: lip n (x) = 0 n{x) e ~ u d^) w h ere 0 is the activity 
parameter, n(x) is the number of objects in the configuration x 
and U s (x) = t. s (xi,Xj) is the interaction potential 
l<i<j<n(x) 
of objects in x. The interaction function between object pairs 
ts(xi,Xj) can be written as follows: 
0 if.Ti~.Tj, 
+00 otherwise. 
(7) 
By assigning an infinite energy to object pair that does not interact 
according to the associated probability is practically zero. 
PLw{9,x,y) = exp j! \e(u:x,y)A(du) 
J Л в (тi\x,y) 
(11) 
where \e is the extended Papangelou intensity which can be writ 
ten as follows: 
Ao(u\x,y) = 0 exp -7dU d {u) - ^ t a (u,Xi) 
\ Xj(zX/x-i^u J 
(12) 
In practice, the simulation step is performed using a Reversible 
Jump Metropolis-Hastings algorithm (Green, 1995) or a multiple 
births and deaths algorithm. Afterwards, the maximization step is 
performed using the Nelder-Mead Simplex algorithm. The SEM 
algorithm is supposed to converge when the parameter vector re 
mains stable. Once the parameters are estimated, we determine 
the most likely object configuration that maximizes the energy 
model thanks to a simulated annealing algorithm.