# A Mutual Reference Shape for Segmentation Fusion and Evaluation

S. Jehan-Besson (a), R. Clouard (b), C. Tilmant (c), A. de Cesare (d), A.
Lalande (e), J. Lebenberg (f), P. Clarysse (g), L. Sarry (h), F. Frouin (i),
M. Garreau (j)  CNRS now with LITO U1288 INSERM Institut Curie Paris,
<EMAIL_ADDRESS>Normandie Univ, UNICAEN, ENSICAEN, CNRS,
GREYC, 14000 Caen, France  Institut Pascal, UMR 6602 UBP CNRS, Clermont-
Ferrand, France  Sorbonne Universités, UPMC Univ. Paris 06, CNRS, INSERM,
Laboratoire d’Imagerie Biomédicale (LIB), 75013, Paris, France  ImViA EA 7535
laboratoire, Université de Bourgogne, France  CEA NeuroSpin UNATI, France Univ
Lyon, INSA-Lyon, Université Lyon 1, UJM-Saint Etienne, CNRS, Inserm, CREATIS
UMR 5220, U1206, F-69621, LYON, France.  Clermont Auvergne Université, CNRS,
Institut Pascal, Clermont-Ferrand, France LITO INSERM U1288 Institut Curie,
Paris, France INSERM U1099, Université de Rennes 1, LTSI, Rennes, F-35000
France

(Written on March 2017)

###### Abstract

This paper proposes the estimation of a mutual shape from a set of different
segmentation results using both active contours and information theory. The
mutual shape is here defined as a consensus shape estimated from a set of
different segmentations of the same object. In an original manner, such a
shape is defined as the minimum of a criterion that benefits from both the
mutual information and the joint entropy of the input segmentations. This
energy criterion is justified using similarities between information theory
quantities and area measures, and presented in a continuous variational
framework. In order to solve this shape optimization problem, shape
derivatives are computed for each term of the criterion and interpreted as an
evolution equation of an active contour. A mutual shape is then estimated
together with the sensitivity and specificity of each segmentation. Some
synthetic examples allow us to cast the light on the difference between the
mutual shape and an average shape. The applicability of our framework has also
been tested for segmentation evaluation and fusion of different types of real
images (natural color images, old manuscripts, medical images).

###### keywords:

Mutual shape , segmentation , variational approaches , segmentation fusion ,
segmentation evaluation , active contours , shape gradients , shape
optimization , average shape , information theory

††journal: ArXiv

##  1 Introduction

Constructing a “consensus” shape from a set of different segmentation results
is an important point when dealing with image segmentation evaluation when no
expert reference is available (evaluation without gold standard). It is also a
key point for an appropriate fusion of several segmentation results in a
single shape. Such a shape must ideally take advantage of the information
provided by each input shape while being robust to outliers. We propose to
tackle the estimation of such a reference shape using information theory
(mutual information and joint entropy) through the definition of a shape
optimization problem. The consensus shape will then be defined as the minimum
of an original criterion based on information theory and area measures and
computed within the framework of active contours and shape gradients. This
shape is called “a mutual shape” and its applicability is tested for image
segmentation evaluation and fusion.

In the context of segmentation evaluation, the estimation of such a reference
shape may be important especially when dealing with large databases of medical
images when the manual delineation of all the frames by an expert becomes a
tedious and time consuming task. The obtained contour is also subject to
inter- and intra-variability and being expert-dependent, it can then not be
considered as an absolute reference. As far as the evaluation without gold
standard is concerned, the STAPLE algorithm (Simultaneous Truth and
Performance Level Estimation) proposed by Warfield et al. [41] is now
classically used in this difficult context. Their algorithm consists in one
instance of the EM (Expectation Maximization) algorithm where the true
segmentation is estimated by maximizing the likelihood of the complete data.
Their pixel-wise approach leads to the estimation of a reference shape
simultaneously with the sensitivity and specificity of each input
segmentation. From these measures, the performance level of each input
segmentation can be estimated and a classification of all the segmentation
entries can be performed. The most recent MAP-STAPLE approach [8] is semi-
local and takes benefit of a small window or patch around the pixel. In this
paper, we rather propose to estimate the reference shape within a continuous
optimization setting by considering such a shape estimation under the umbrella
of shape optimization tools [14] and deformable models [21]. Indeed, the
computation of a reference shape can be advantageously modeled as the optimum
of a well chosen energy criterion and estimated by a shape gradient descent
that corresponds to the deformation of an active shape. Moreover, we propose a
new theoretical criterion based on information theory that allows to well
understand the behaviour of our reference shape.

Let us also note that shape optimization algorithms have already been proposed
in order to compute shape averages [7, 43] or more recently median shapes [2]
by minimizing different shape metrics like the Hausdorff distance in [7] or
the symmetric area difference between shapes in [43]. Some other approaches
also take advantage of well-appropriate distances between level-set shapes
(see for example [28]). Comparing with these previous approaches, our goal is
quite different since our aim is to compute a consensus shape from $N$ input
segmentations.

The main contribution of this paper is then to propose a new theoretical model
to carry out the estimation of a consensus or reference shape from several
segmentation entries using active contours and shape gradients. This
theoretical model is based on information theory and justified using the
analogies between information theory and area measures. In order to estimate
what we call a “mutual shape”, we then propose to maximize the mutual
information between the $N$ input segmentations while minimizing the joint
entropy. Such a statistical criterion can be interpreted as a robust measure
of the symmetric area difference. The minimization is performed through the
computation of the evolution equation of an active contour. This evolution
equation is computed using advanced shape derivation tools. In order to
perform such a derivation, the criterion must be expressed in a continuous
settings and non parametric probability density functions are estimated using
Kernel methods [15]. In this variational setting, we also propose to add a
classical regularization term based on the curvature of the deformable
contour. Such a term is weighted using a regularization parameter that
controls the smoothness of the obtained contour. The advantage of this
formalism is to make explicitly appear, in the criterion to minimize, the
domain and the associated contour. The criterion is then easier to understand
and interpret and some geometrical and photometric priors could be directly
added in the criterion to minimize. The derivation is directly performed
according to the domain using shape derivation tools.

The proposed algorithm is first experimented on a synthetic example that
allows to understand the differences between a classic average shape based on
a symmetric area minimization [43], a simple majority voting shape and the
proposed mutual shape. It is also evaluated for segmentation fusion and
evaluation on different images : color real natural images, old manuscripts or
medical images, in order to show the genericity of this framework. The first
application concerns segmentation evaluation and fusion on a real color
natural image using segmentations from the Berkeley database [29]. The second
application is dedicated to text segmentation in old manuscripts and we
propose two main examples of segmentation. One of them takes benefit of the
DIBCO database [33]. The last application is devoted to segmentation fusion
and evaluation of different delineation methods of the left ventricular cavity
in Magnetic Resonance Imaging (MRI). For this application, we propose to
compare the mutual shape to the reference algorithm STAPLE [41] classically
used for segmentation fusion and evaluation in medical images.

In section 2, we present the problem statement and in section 3, the proposed
criterion for the estimation of the mutual shape. The criterion is then
estimated in a continuous framework and expressed using domain or contour
integrals in section 4. Such a continuous criterion can then be derived using
shape optimization tools in order to compute the mutual shape (see section 5).
Experimental results on synthetic examples are detailed in section 6 and the
different applications in section 7.

##  2 Problem statement

Let $\mathcal{U}$ be a class of domains (open regular bounded sets, i.e.
$\mathcal{C}^{2}$) of $\mathbb{R}^{d}$ (with $d=2$ or $3$). In this paper
theoretical results are stated for $d=2$ or $d=3$ but the experimental results
are conducted on 2D-images. We denote by $\Omega_{i}$ an element of
$\mathcal{U}$ of boundary $\partial\Omega_{i}$. We consider
$\\{\Omega_{1},...,\Omega_{n}\\}$ a family of $n$ shapes where each shape
corresponds to the segmentation of the same unknown object $O$ in a given
image. The image domain is denoted by $\Omega\in\mathbb{R}^{d}$. Our aim is to
compute a reference shape $\mu$ that can closely represent the true object $O$
(Fig.1). We propose to define the problem through a statistical representation
of shapes embedded in an information theory criterion. Let us first recall the
main shape representation models and criteria proposed in the literature.

Figure 1: Diagram of the problem statement: evaluation of a reference shape
$\mu$ from a set of $n$ segmented shapes of the same object.

###  2.1 Shape representation

The computation of a reference shape is closely linked to the choice of a
representation. An analytical representation may be used as in [22] where the
authors propose a statistical study of shapes by representing them as a finite
number of points. Some authors prefer to choose an implicit representation of
shapes which avoids the parametrization step. For example, in [2, 6] shapes
are represented using their characteristic function as follows:

$\displaystyle d_{i}({\mathbf{x}})=\left\\{\begin{array}[]{c c
c}1&\text{if}&{\mathbf{x}}\in\Omega_{i}\\\
0&\text{if}&{\mathbf{x}}\,\not\in\Omega_{i}\end{array}\right.$ (3)

where ${\mathbf{x}}\in\Omega$ is the location of the pixel within the image.
We denote by $\overline{\Omega_{i}}$ the complementary shape of $\Omega_{i}$
in $\Omega$ with $\Omega_{i}\cup\overline{\Omega_{i}}=\Omega$. One may also
take advantage of the distance function associated to each shape. In [28] the
authors propose to perform a principal component analysis on shapes in order
to provide a statistical shape prior. In the same vein, some statistical shape
priors have been proposed by [13, 32] using this implicit representation.

More recently shapes have been represented using Legendre moments in order to
define shape priors for segmentation using active contours [17]. This
representation can also be easily included in a variational setting [17, 27].

We may also consider that each shape is a realization of a random variable.
Such a representation has been introduced in [41] in order to evaluate a
reference shape in a statistical framework, in [40] for the morphological
exploration of shape spaces and statistics, and also in [18, 23] for image
segmentation using information theory. In this paper, we take advantage of
this statistical representation that appears to be well adapted to the
definition of a statistical criterion. The shape is represented through a
random variable $D_{i}$ whose observation is the characteristic function
$d_{i}$ defined in (3). The reference shape $\mu$ is also represented through
an unknown random variable $T$ with the associated characteristic function
$t({\mathbf{x}})=1$ if ${\mathbf{x}}\in\mu$ and $t({\mathbf{x}})=0$ if
${\mathbf{x}}\in\overline{\mu}$.

###  2.2 Definition of average shapes

In the literature, average shapes are defined through the minimization of the
sum of the distances of the unknown shape $\mu$ to each shape $\Omega_{i}$ as
follows:

$\displaystyle\mu=\arg\min_{\mu^{*}}\sum_{i=1}^{n}d(\Omega_{i},\mu^{*})$ (4)

Of course, the definition of the distance $d$ is crucial and may lead to
different results and average shapes. For example, an average shape can be
computed by minimizing the area of the symmetric differences [43] using
$d(\Omega_{i},\mu):=|\Omega_{i}\triangle\mu|$ where $|.|$ stands for the
cardinal of the considered domain. In a continuous optimization framework, the
criterion to minimize can be expressed as follows:

$\displaystyle
SD(\mu)=\sum_{i=1}^{n}|\Omega_{i}\triangle\mu|=\sum_{i=1}^{n}\left(\int_{\mu}(1-d_{i}({\mathbf{x}}))d{\mathbf{x}}+\int_{\overline{\mu}}d_{i}({\mathbf{x}})d{\mathbf{x}}\right)$
(5)

In [6, 7], the authors prefer to introduce the Hausdorff distance to perform
shape warping while in [2], the authors modify the previous criterion in order
to compute a median shape. In addition to the previous works, we can also cite
[40] where the authors propose to explore shape spaces using mathematical
morphology. The optimal shape is computed using a watershed performed on the
squared sum of the distance functions or using a morphological computation of
a median set. Another class of algorithms was proposed for the estimation of
an unknown shape from multiple channels (color or multimodal segmentation). We
can cite the work of Chan et al. [5] or the multimodal segmentation approaches
proposed in [18, 23]. These works were not designed at first for segmentation
evaluation or fusion but they are worth mentioning because they propose to
treat the different channels in a single criterion (may also be useful for
information fusion). Moreover in [18, 23], some information theory quantities
are used. Our work is different especially due to the fact that we consider
both the maximization of mutual information coupled with the minimization of
joint entropies and the joint estimation of evaluation quantities (sensitivity
and specificity measures).

##  3 Proposition of a criterion for the estimation of a mutual shape

Our goal is here to mutualize the information given by each segmentation to
define a consensus or reference shape. Such a shape cannot be considered as a
simple average shape.

In this context, we propose to take advantage of the analogies between
information measures (mutual information, joint entropy) and area measures. As
previously mentioned, $D_{i}$ represents the random variable associated with
the characteristic function $d_{i}$ of the shape $\Omega_{i}$ and $T$ the
random variable associated with the characteristic function $t$ of the
reference shape $\mu$. Using these notations, $H(D_{i},T)$ represents the
joint entropy between the variables $D_{i}$ and $T$, and $I(D_{i}/T)$ the
mutual information. We then propose to minimize the following criterion :

$\displaystyle
E(T)=\sum_{i=1}^{n}\left(H(D_{i},T)-I(D_{i},T)\right)=JH(T)+MI(T)$ (6)

where the sum of joint entropies is denoted by
$JH(T)=\sum_{i=1}^{n}H(D_{i},T)$ and the sum of mutual information by
$MI(T)=-\sum_{i=1}^{n}I(D_{i},T)$.

Introducing this criterion can be justified by the fact that
$\varphi(D_{i},T)=(H(D_{i},T)-I(D_{i},T))$ is a metric which satisfies the
following properties :

  1. 1.

$\varphi(X,Y)\geq 0$

  2. 2.

$\varphi(X,Y)=\varphi(Y,X)$

  3. 3.

$\varphi(X,Y)=0$ if and only if $X=Y$

  4. 4.

$\varphi(X,Y)+\varphi(Y,Z)\geq\varphi(X,Z)$

Indeed, we can show easily that $\varphi(D_{i},T)=H(T/D_{i})+H(D_{i}/T)$ using
the following classical relations between the joint entropy and the
conditional entropy and the mutual information and the conditional entropy :

$\displaystyle H(D_{i},T)=H(D_{i})+H(T/D_{i})$ $\displaystyle
I(D_{i},T)=H(D_{i})-H(D_{i}/T)$

Finally, $H(T/D_{i})+H(D_{i}/T)$ is shown to be a metric that satisfies the
four properties above [12]. We then minimize a sum of distances between
$D_{i}$ and $T$ expressed using information theory quantities.

Moreover, we can give a further interesting geometrical interpretation of the
proposed criterion. We propose to take advantage of the analogies between
information measures (mutual information, joint entropy) and area measures. In
[35, 42], it is shown that Shannon’s information measures can be interpreted
in terms of area measures as follows:

$\displaystyle
H(D_{i},T)=\operatorname{mes}(\tilde{D_{i}}\cup\tilde{T})\quad\text{and}\quad
I(D_{i},T)=\operatorname{mes}(\tilde{D_{i}}\cap\tilde{T}),$ (7)

with $\tilde{X}$ the abstract set associated with the random variable $X$ and
the term $\operatorname{mes}$ corresponds to a signed measure defined on an
algebra of sets with values in $]-\infty,+\infty[$. The signed measure must
satisfy $\operatorname{mes}(\emptyset)=0$ and
$\operatorname{mes}(\bigcup_{k=1}^{n}A_{k})=\sum_{k=1}^{n}\operatorname{mes}(A_{k})$
for any sequence $\\{A_{k}\\}_{k=1}^{n}$ of disjoint sets. Each quantity can
then be viewed as an operation on the sets (Fig.2). These properties help us
to better understand the role of each term in the criterion to optimize.

Figure 2: Mutual information and joint entropy as area measures

Indeed, when estimating a classic average shape using the criterion (5), one
performs the minimization of the sum of the union of the shapes $\Omega_{i}$
with $\mu$ while maximizing the sum of the intersection between the same
shapes. By analogy with this criterion, we minimize a measure of the union
while maximizing a measure of the intersection through the use of information
quantities. In other words, the sum of the joint entropies (union of sets)
will be minimized while the sum of the mutual information quantities
(intersection) will be maximized. The proposed criterion can then be
interpreted as a statistical measure of the area of the symmetric difference
which is really interesting for the estimation of a consensus shape.

##  4 Expression of the criterion in a continuous framework

In order to take advantage of the previous statistical criterion (6) within a
continuous shape optimization framework, we propose to express the joint and
conditional probability density functions according to the reference shape
$\mu$. This step is detailed in this section for both the mutual information
and the joint entropy.

###  4.1 Maximization of mutual information (MI)

Here we try to express $MI(T)=-\sum_{i=1}^{n}I(D_{i},T)$ in a continuous
setting according to the unknown shape $\mu$. In order to simplify the
criterion, we use the classic relation between mutual information and
conditional entropy: $I(D_{i},T)=H(D_{i})-H(D_{i}/T)$. Since $H(D_{i}/T)\geq
0$ and $H(D_{i})$ is independent of $T$, we will rather minimize
$\sum_{i=1}^{n}H(D_{i}/T)$. Denoting by $t$ and $d_{i}$ the observations of
the random variables $T$ and $D_{i}$, the conditional entropy of $D_{i}$
knowing $T$ can be written as follows:

$\displaystyle
H(D_{i}/T)=-\sum_{t\in\\{0,1\\}}\left[p(t)\sum_{d_{i}\in\\{0,1\\}}p(d_{i}/t)\log(p(d_{i}/t))\right],$
(8)

with $p(T=t)=p(t)$ and $p(D_{i}=d_{i}/T=t)=p(d_{i}/t)$.
The conditional probability $p(d_{i}=1/t=1)$ corresponds to the sensitivity
parameter $p_{i}$ (true positive fraction):

$\displaystyle
p_{i}(\mu)=p(d_{i}=1/t=1)=\frac{1}{|\mu|}\int_{\mu}K(d_{i}({\mathbf{x}})-1)d{\mathbf{x}}.$
(9)

where the function $K$ represents a Gaussian Kernel of $0$-mean and variance
$\sigma$. This function allows a rigorous application of the shape derivation
tools due to the fact that the function under the integral is differentiable.
In this paper, we choose a very small $\sigma=0.1$ for all the experimental
studies.
The conditional probability $p(d_{i}=0/t=0)$ corresponds to the specificity
parameter $q_{i}$ (true negative fraction) :

$\displaystyle
q_{i}(\mu)=p(d_{i}=0/t=0)=\frac{1}{|\overline{\mu}|}\int_{\overline{\mu}}K(d_{i}({\mathbf{x}}))d{\mathbf{x}}.$
(10)

In the rest of the paper, for the sake of simplicity, $p_{i}(\mu)$ is replaced
by $p_{i}$ and $q_{i}(\mu)$ by $q_{i}$. The random variable $T$ takes the
value $1$ with a probability $p(t=1)=|\mu|/|\Omega|$ and $0$ with a
probability $p(t=0)=|\overline{\mu}|/|\Omega|$. The $MI$ criterion can then be
expressed according to $\mu$:

$\displaystyle MI(\mu)=$ $\displaystyle-$
$\displaystyle\sum_{i=1}^{n}\Big{[}\frac{|\mu|}{|\Omega|}\left((1-p_{i})\log\left(1-p_{i}\right)+p_{i}\log
p_{i}\right)$ (11) $\displaystyle+$
$\displaystyle\frac{|\overline{\mu}|}{|\Omega|}\left(q_{i}\log
q_{i}+\left(1-q_{i}\right)\log\left(1-q_{i}\right)\right)\Big{]}$

The parameters $p_{i}$ and $q_{i}$ depend explicitly on $\mu$, which must be
taken into account in the optimization process. Indeed if $\mu$ is updated in
an iterative process, the parameters $p_{i}$ and $q_{i}$ must also be updated
which implies a joint estimation of these quantities with the unknown mutual
shape.

###  4.2 Minimization of joint entropy

Let us now express, according to $\mu$ and in a continuous setting, the sum of
the joint entropies $JH(T)=\sum_{i=1}^{n}H(D_{i},T)$. The following expression
of the joint entropy is considered:

$\displaystyle
H(D_{i},T)=-\sum_{t\in\\{0,1\\}}\sum_{d_{i}\in\\{0,1\\}}p(d_{i},t)\log\left(p(d_{i},t)\right),$
(12)

with $p(D_{i}=d_{i},T=t)=p(d_{i},t)$.
The following estimates for the joint probabilities are then used ($a=0$ or
$a=1$):

$\displaystyle
p(d_{i}=a,t=1)=\frac{1}{|\Omega|}\int_{\mu}\left(K(d_{i}({\mathbf{x}})-a)\right)d{\mathbf{x}},$
$\displaystyle
p(d_{i}=a,t=0)=\frac{1}{|\Omega|}\int_{\overline{\mu}}\left(K(d_{i}({\mathbf{x}})-a)\right)d{\mathbf{x}}.$
(13)

where the function $K$ represents a Gaussian Kernel of $0$-mean and variance
$\sigma$. The criterion to minimize is now denoted by $JH(\mu)$ and can be
written as follows:

$\displaystyle JH(\mu)$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{n}\Big{[}p(d_{i}=1,t=1)\log(p(d_{i}=1,t=1))$
$\displaystyle-$ $\displaystyle p(d_{i}=1,t=0)\log(p(d_{i}=1,t=0))$
$\displaystyle-$ $\displaystyle p(d_{i}=0,t=1)\log(p(d_{i}=0,t=1))$
$\displaystyle-$ $\displaystyle p(d_{i}=0,t=0)\log(p(d_{i}=0,t=0))\Big{]}$

where $p(d_{i}=a,t=1)$ and $p(d_{i}=a,t=0)$ depends on $\mu$ as expressed in
equations (13).

###  4.3 Continuous expression of the criterion

Using the two previous sections, we can express the global criterion to
minimize according to $\mu$ as follows:

$\displaystyle E(\mu)$ $\displaystyle=$ $\displaystyle JH(\mu)+MI(\mu)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\Big{[}-p(d_{i}=1,t=1)\log(p(d_{i}=1,t=1))$
$\displaystyle-$ $\displaystyle p(d_{i}=1,t=0)\log(p(d_{i}=1,t=0))$
$\displaystyle-$ $\displaystyle p(d_{i}=0,t=1)\log(p(d_{i}=0,t=1))$
$\displaystyle-$ $\displaystyle p(d_{i}=0,t=0)\log(p(d_{i}=0,t=0))$
$\displaystyle-$
$\displaystyle\sum_{i=1}^{n}\Big{[}\frac{|\mu|}{|\Omega|}\left((1-p_{i})\log\left(1-p_{i}\right)+p_{i}\log
p_{i}\right)$ $\displaystyle+$
$\displaystyle\frac{|\overline{\mu}|}{|\Omega|}\left(q_{i}\log
q_{i}+\left(1-q_{i}\right)\log\left(1-q_{i}\right)\right)\Big{]}$

In this given form, the minimization of such a criterion can be considered
using active contours and shape gradients as detailed in the following
section.

##  5 Optimization using shape gradients

In order to compute a local minimum of the criterion $E$ defined in (4.3), we
propose to take advantage of the framework developed in [1] which is based on
the shape optimization tools proposed in [14, Chap.8]. The main idea is to
deform an initial curve (or surface) towards the boundaries of the region of
interest.

Formally, the contour then evolves according to the following Partial
Differential Equation (PDE):

$\displaystyle\frac{\partial\Gamma(z,\tau)}{\partial\tau}=v({\mathbf{x}},\mu){\mathbf{N}}({\mathbf{x}})$
(16)

where $\Gamma(z,\tau)$ is the evolving curve, $z$ a parameter of the curve,
$\tau$ the evolution parameter, $v({\mathbf{x}},\mu)$ the amplitude of the
velocity in ${\mathbf{x}}=\Gamma(z,\tau)$ directed along the normal of the
curve ${\mathbf{N}}({\mathbf{x}})$. The evolution equation and more
particularly the velocity $v$ must be computed in order to make the contour
evolve towards an optimum of the energy criterion. From an initial curve
$\Gamma_{0}$ defined by the user, we will have
$\lim\limits_{\tau\rightarrow\infty}\Gamma(\tau)=\mu$ at convergence of the
process.

The main issue lies in the computation of the velocity $v$ in order to find
the unknown shape $\mu$ at convergence. This term is deduced from the
derivative of the criterion according to the shape. The method of derivation
is explained in details in [1] and is based on shape derivation principles
developed formally in [14, 38]. For completeness, we recall some useful
definitions and theorems.

###  5.1 Main mathematical tools

The following theorem is the central theorem for derivation of integral
domains of the form $\int_{\mu}k({\mathbf{x}},\mu)\,d{\mathbf{x}}$. It gives a
general relation between the Eulerian derivative and the shape derivative for
region-based terms.

###### Theorem 1

Let $\Omega$ be a $C^{1}$ domain in $\mathbb{R}^{n}$ and ${\mathbf{V}}$ a
$C^{1}$ vector field. Let $k$ be a $C^{1}$ function. The functional
$J(\mu)=\int_{\mu}k({\mathbf{x}},\mu)\,d{\mathbf{x}}$ is differentiable and
its Eulerian derivative in the direction of ${\mathbf{V}}$ is the following:

$\displaystyle<J^{\prime}(\mu),{\mathbf{V}}>=\int_{\mu}k_{s}({\mathbf{x}},\mu)\,d{\mathbf{x}}-\int_{\partial\mu}k({\mathbf{x}},\mu)({\mathbf{V}}\cdot{\mathbf{N}})\,d{\mathbf{a}}$
(17)

where $k_{s}$ is the shape derivative of $k$ defined by
$k_{s}({\mathbf{x}},\mu)=\lim_{\tau\rightarrow
0}\frac{k({\mathbf{x}},\mu(\tau))-k({\mathbf{x}},\mu)}{\tau}$. The term
${\mathbf{N}}$ denotes the unit inward normal to $\partial\mu$ and
$d{\mathbf{a}}$ its area element (in $\mathbb{R}^{2}$, we have
$d{\mathbf{a}}=ds$ where $s$ stands for the arc length).

The Eulerian derivative of $J$ in the direction ${\mathbf{V}}$ is defined as

$<J^{\prime}(\mu),{\mathbf{V}}>=\lim_{\tau\rightarrow
0}\frac{J(\mu(\tau))-J(\mu)}{\tau}$

if the limit exists, with $\mu(\tau)=T_{\tau}({\mathbf{V}})(\mu)$ the
transformation of $\mu$ through the vector field ${\mathbf{V}}$. The proof of
the theorem can be found in [14].

###  5.2 Methodology for the computation of the evolution equation

The following proposition gives us a way to compute the evolution equation of
the active contour when the Eulerian derivative can be expressed as an
integral over the boundary of the domain.

###### Proposition 1

Let us consider that the shape derivative of the criterion $J(\mu)$ in the
direction ${\mathbf{V}}$ may be written in the following way:

$<J^{\prime}(\mu),{\mathbf{V}}>=-\int_{\partial\mu}v({\mathbf{x}},\mu)({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$
(18)

Interpreting this equation as the $L^{2}$ inner product on the space of
velocities, the straightforward choice in order to minimize $J(\mu)$ consists
in choosing ${\mathbf{V}}=v{\mathbf{N}}$ for the deformation. We can then
deduce that, from an initial contour $\Gamma_{0}$, the boundary $\partial\mu$
can be found at convergence of the following evolution equation:

$\displaystyle\frac{\partial\Gamma}{\partial\tau}=v({\mathbf{x}},\mu)\,{\mathbf{N}}$
(19)

where $v$ is the velocity of the curve and $\tau$ the evolution parameter.

The shape derivatives of the criteria $SD(\mu)$ (5), $MI(\mu)$ (11) and
$JH(\mu)$ (4.2), can be written in the form (18) which allows us to find some
geometrical PDEs of the form (19) for each criterion. The derivation is
developed thereafter.

###  5.3 Shape derivatives

This paragraph details the shape derivatives of the criteria $SD(\mu)$ (5),
$MI(\mu)$ (11) and $JH(\mu)$ (4.2). Proofs of the two main new theorems 3 and
4 are given in the appendix of the paper.

###### Theorem 2

The shape derivative in the direction ${\mathbf{V}}$ of the functional
$SD(\mu)$ given in (5) is:

$<SD^{\prime}(\mu),{\mathbf{V}}>=-\int\limits_{\Gamma}v_{SD}\,({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$

with the velocity :

$\displaystyle\quad v_{SD}=\sum_{i=1}^{n}(1-2\,d_{i}({\mathbf{x}})).$ (20)

The computation of the shape derivative of $MI(\mu)$ is more complex because
the functions inside the integrals depend on $\mu$.

###### Theorem 3

The shape derivative in the direction ${\mathbf{V}}$ of the functional
$MI(\mu)$ defined in (11) is:

$<MI^{\prime}(\mu),{\mathbf{V}}>=-\int_{\Gamma}v_{MI}\,({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$

with the velocity

$\displaystyle v_{MI}=$
$\displaystyle\frac{1}{|\Omega|}\sum_{i=1}^{n}\Big{[}(p_{i}-K(d_{i}-1))\log\left(\frac{pi}{1-p_{i}}\right)$
$\displaystyle-$
$\displaystyle(q_{i}-K(d_{i}))\log\left(\frac{qi}{1-q_{i}}\right)$
$\displaystyle+$ $\displaystyle q_{i}\log q_{i}+(1-q_{i})\log(1-q_{i})$
$\displaystyle+$ $\displaystyle p_{i}\log p_{i}+(1-p_{i})\log(1-p_{i})\Big{]}$

The computation of the shape derivative of $JH(\mu)$ is also complex and leads
to the following theorem:

###### Theorem 4

The shape derivative in the direction ${\mathbf{V}}$ of the functional
$JH(\mu)$ defined in (4.2) is:

$<JH^{\prime}(\mu),{\mathbf{V}}>=-\int_{\Gamma}v_{JH}\,({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$

The velocity $v_{JH}$ is given by the following equation:

$\displaystyle v_{JH}$ $\displaystyle=$
$\displaystyle\frac{-1}{|\Omega|}\sum_{i=1}^{n}\Big{[}K(d_{i}-1)\log\left(\frac{p(d_{i}=1,t=1)}{p(di=1,t=0)}\right)$
(22) $\displaystyle+$ $\displaystyle
K(d_{i})\log\left(\frac{p(d_{i}=0,t=1)}{p(di=0,t=0)}\right)\Big{]}.$

where $v_{JH}$ is directed along ${\mathbf{N}}$.

###  5.4 Global evolution equations for the different criteria

A standard regularization term is added in the criterion to minimize in order
to favor smooth shapes. This term corresponds to the minimization of the curve
length and is defined by $Reg(\mu)=\int_{\partial\mu}ds$. It is balanced with
a positive coefficient $\lambda$ in the criterion and leads to the following
velocity in the evolution equation:

$\displaystyle v_{Reg}=\kappa$ (23)

where $\kappa$ is the curvature of the contour $\Gamma(\tau)$.

Finally, we propose to define our mutual reference shape through the
minimization of a global criterion called $J_{IT}$ (Information Theoretic
criterion):

$\displaystyle J_{IT}(\mu)=JH(\mu)+MI(\mu)+\lambda Reg(\mu).$ (24)

In order to minimize this criterion, the following evolution equation is used:

$\displaystyle\left(\frac{\partial\Gamma}{\partial\tau}\right)_{IT}=\left(v_{JH}+v_{MI}+\lambda
v_{Reg}\right)\bf{N}$ (25)

where $v_{MI}$, $v_{JH}$ and $v_{Reg}$ are defined respectively in equations
(3), (22) and (23). The term $\bf{N}$ designates the inward unit normal of the
active contour. In the experimental results, the mutual reference shape is
also compared to the average shape (SD) that corresponds to the minimization
of the following criterion:

$\displaystyle J_{SD}(\mu)=SD(\mu)+\lambda Reg(\mu).$ (26)

In order to minimize this criterion, the following evolution equation is
applied:

$\displaystyle\left(\frac{\partial\Gamma}{\partial\tau}\right)_{SD}=\left(v_{SD}+\lambda
v_{Reg}\right)\bf{N}$ (27)

where $v_{SD}$ and $v_{Reg}$ are defined respectively in equations (20) and
(23). These velocities are directed along the unit inward normal $\bf{N}$ of
the active contour.

Note also that using this formalism, some other prior information (photometric
or geometric) can be inserted by adding some additional velocities in the PDE.
For example, we may take advantage of the tools developed in [4, 17, 27, 32].

###  5.5 Implementation of the active contour

As far as the numerical implementation is concerned, we use the level set
method [30]. The key idea is to introduce an auxiliary function
$U({\mathbf{x}},\tau)$ such that $\Gamma(\tau)$ is the zero level set of $U$.
The function $U$ is often chosen to be the signed distance function of
$\Gamma(\tau)$. The evolution equation then becomes:

$\displaystyle\frac{\partial U}{\partial\tau}=F|\nabla U|.$ (28)

The velocity is chosen as $F=v_{JH}+v_{MI}+\lambda v_{Reg}$ for the estimation
of the mutual shape and $F=v_{SD}+\lambda v_{Reg}$ for the estimation of the
SD shape. This method is accurate and allows to automatically handle the
topological changes of the initial curve. However, the same evolution
equations could be implemented using faster implementation algorithms such as
B-splines [34]. Convex optimization methods [3] may perhaps be interesting but
the criterion is not convex and some assumptions are needed before a direct
application of these methods.

##  6 Experimental results on a synthetic example

The behavior of our mutual shape estimation is first tested on a synthetic
example. The mutual shape, the classic average shape and a simple majority
voting approach are compared. We also study the joint evolution of the
sensitivity and specificity parameters.

###  6.1 Difference between a mutual shape and a classic average shape

In this section, a test sequence consisting of different segmentations of a
lozenge (Fig.3) was built. The first entry is the true segmentation mask, the
other entries represent the segmentation of $1/4$ of the true lozenge
(Fig.3.b).

|
---|---
(a) | (b)
Figure 3: The image to segment is given in (a) and the different segmentation
entries (masks) for this image are given in (b).

When computing the average of the different characteristic functions using the
formula : $\sum_{i=1}^{n}d_{i}/n$, we remark (Fig.4.b) that some masks share
an intersection. Indeed the values of the average image belong to the interval
$[0,0.6]$. The value $0$ corresponds to black points in Fig.4.a and the value
$0.6$ corresponds to white points. We then binarize this average image $I_{A}$
in an image named $I_{AT}$ displayed in (Fig.4.b). If $I_{A}({\mathbf{x}})\geq
0.5$ then $I_{AT}=0$ (black points) and if $I_{A}({\mathbf{x}})<0.5$ then
$I_{AT}=255$ (white points). This procedure gives us a simple majority voting
procedure. The result is the black line inside the lozenge.

|
---|---
(a) | (b)
Figure 4: The average image $I_{A}$ (a) and the corresponding binarized
average image $I_{AT}$ (b) of the masks of the Fig.3(b) (simple majority
voting procedure).

We then use the evolution equations of both the mutual shape (25) and of the
SD shape (27). The initial contour is chosen as a circle including the lozenge
(Fig.5.a and Fig.6.a). The mutual shape algorithm is able to recover the whole
lozenge and is then different from a classic average shape (see Fig.5 and
Fig.6). The curve evolves and segments the whole lozenge by an iterative
process (images resulted from different iterations in Fig.5.b and Fig.5.c).
The final contour is given in Fig.5.d. The mutual shape is compared to a shape
average computed using the minimization of the classic symmetrical difference
(criterion $J_{SD}$ with evolution equation (27)). The evolution is given in
Fig.6. In this case, the final contour is similar to the result obtained by
computing a binarized mean $I_{AT}$ (Fig.4.b) since it corresponds to a line
due to the small overlap between masks 2 and 5 (Fig.3.b). The same small value
is taken for the regularization parameter $\lambda$ in order to give a higher
importance to the data term.

|  |  |
---|---|---|---
(a) Initial contour | (b) It. 80 | (c) It. 140 | (d) Mutual shape
Figure 5: Evolution using the mutual shape (evolution equation (25) with $\lambda=10$). In the first image (a), the initial contour is in white (circle) and the other white lines represent the boundaries of the different segmentation entries. Intermediate results obtained from $80$ and $140$ iterations are displayed in images (b) and (c) and the final estimated mutual shape in (d) (240 iterations). |  |  |
---|---|---|---
(a) Initial contour | (b) It. 300 | (c) It. 400 | (d) SD shape
Figure 6: Evolution using the SD shape (evolution equation (27) with
$\lambda=10$). In the first image, the initial contour is in white (circle)
and the other white lines represent the boundaries of the different
segmentation entries. Intermediate results obtained from $300$ and $400$
iterations are displayed in images (b) and (c) and the final estimated SD
shape in (d) (600 iterations).

###  6.2 Difference between the mutual shape and the union of the masks

An outlier (Fig.7.a) was introduced in the initial sequence of masks in order
to test the robustness of the mutual shape estimation. Indeed, our goal is to
test that the mutual shape is also different to a simple union of the
different masks. In Fig.7, the different steps of the evolution of the contour
are displayed. The final contour (Fig.7.d) fits the lozenge and excludes the
outlier from the final contour.

|  |  |  |
---|---|---|---|---
(a) Input outlier | (b) Initial | (c) It. 100 | (d) It. 380 | (e) Mutual shape
Figure 7: Introduction of an outlier (a) in the initial sequence of masks
(Fig.3.b) and estimation of the mutual shape (evolution equation (25) with
$\lambda=10$). In the image (b), the initial contour is in white and the other
white contours and lines represent the different boundaries of the initial
masks (the segmentation entries and the outlier).

###  6.3 Joint evolution of the sensitivity and specificity parameters

When the active contour evolves using the evolution equation (25), the
parameters $p_{i}$ and $q_{i}$ are estimated jointly with the mutual shape as
proposed in STAPLE [41]. The different values of these parameters along the
evolution of the curve are given in Table 1. These results are obtained using
masks displayed in the first row of this Table. According to the final values
reported in Table 1, we can conclude that the best segmentation corresponds to
the shape $1$ with $p_{1}=1$ and $q_{1}=1$ and that the shape $6$ is an
outlier since the sensitivity coefficient is equal to $0$. The other
segmentations correspond to one quarter of the lozenge which leads to a
sensitivity parameter around the value of $0.25$. Note that the initial values
of $p_{i}$ and $q_{i}$ are computed directly using the initial contour.

We can notice that the specificity parameter $q_{i}$ is less relevant. Indeed
this parameter is estimated using the external domain ($\bar{\mu}$) and is
then estimated using a higher number of pixels. It should be normalized in
order to be comparable to the $p_{i}$ value. One solution consists in the
selection of a smaller working area (a mask that includes the union of masks
chosen in order to get two regions with a comparable size).

|  |  |  |  |  |
---|---|---|---|---|---|---
Iterations | mask 1 | mask 2 | mask 3 | mask 4 | mask 5 | mask 6
It. 0 | $p_{1}=0.35$ | $p_{2}=0.09$ | $p_{3}=0.09$ | $p_{4}=0.08$ | $p_{5}=0.09$ | $p_{6}=0.15$
(Fig. 7.b) | $q_{1}=1$ | $q_{2}=1$ | $q_{3}=1$ | $q_{4}=1$ | $q_{5}=1$ | $q_{6}=1$
It. 100 | $p_{1}=0.60$ | $p_{2}=0.15$ | $p_{3}=0.15$ | $p_{4}=0.13$ | $p_{5}=0.16$ | $p_{6}=0.27$
(Fig. 7.c) | $q_{1}=1$ | $q_{2}=1$ | $q_{3}=1$ | $q_{4}=1$ | $q_{5}=1$ | $q_{6}=1$
Final | $p_{1}=1$ | $p_{2}=0.24$ | $p_{3}=0.26$ | $p_{4}=0.22$ | $p_{5}=0.27$ | $p_{6}=0$
(Fig. 7.e) | $q_{1}=1$ | $q_{2}=1$ | $q_{3}=1$ | $q_{4}=1$ | $q_{5}=1$ | $q_{6}=0.93$
Table 1: Joint evolution of the contour and of the sensitivity and specificity
parameters $p_{i}$ and $q_{i}$ for the masks 1 to 6. The values correspond to
the evolution of the contour displayed in Fig.7 (initial contour, iteration
100 and final contour).

##  7 Experimental results on real images

In this section, our aim is to provide a variety of examples where the
proposed mutual shape can be valuable. Indeed the theoretical framework
proposed above is generic and can be applied to different images, modalities,
shapes and applications. First of all, in subsection 7.1, we provide a simple
example on a real color image from the Berkeley database [29] to show the
robustness of our estimation to an outlier, the accuracy of the obtained
contour and the relevance of the classification performed using $p_{i}$ and
$q_{i}$. As already mentioned, the implementation is performed using the level
set method which automatically handles topological changes. Therefore, we then
apply the estimation of the mutual shape for complicated shapes composed of
multiple separated components such as the text in old documents. In the
subsection (7.2), we give two examples : the first one is dedicated to the
fusion of very simple binarization techniques while the second one performs
fusion and evaluation of real automatic binarization methods from the DIBCO
challenge [33]. In the subsection 7.3, we propose to test the mutual shape for
the fusion and the evaluation without gold standard of different segmentation
methods or expert delineations of the left ventricle in cardiac magnetic
resonance images (cardiac MRI) and notably expert segmentations. This
estimated mutual shape is compared to the classical STAPLE estimation [41] and
evaluation results are analysed on the basis of some previous works on
evaluation without gold standard [26].

Let us note that the parameter $\lambda$ is chosen small ($1$ to $10$) for non
convex shapes and may be chosen higher for convex shapes ($10$ to $100$). In
this last case, it can help to get a more regularized contour.

###  7.1 Application on a real natural color image

The estimation of such a mutual shape is first tested for the unsupervised
evaluation of segmentation methods of real images. The object of interest is
the tiger of the image displayed in Fig.9. We then extract the object of
interest from the different segmented images proposed in the Berkeley
database[29]. The different segmentation entries $m_{1}$ to $m_{5}$ are given
in Fig.8 and we add an outlier $m_{6}$, which corresponds to the segmentation
of the tree behind the tiger, to the five main segmentation entries.

|  |
---|---|---
$m_{1}$ | $m_{2}$ | $m_{3}$
|  |
$m_{4}$ | $m_{5}$ | $m_{6}$
Figure 8: The different segmentation masks $m_{1}$ to $m_{5}$ and an outlier $m_{6}$ are taken as segmentation entries for the mutual shape estimation.  |
---|---
(a) Initial contour | (b) Iteration 200
|
(c) Estimated mutual shape (final contour) | (d) Final segmentation
Figure 9: Evolution of the active contour (in white, with $\lambda=10$), from the initial contour (a), for the estimation of the mutual shape using the masks $m_{1}$ to $m_{6}$ (Fig.8). The final mutual shape is given in (c) and the segmented tiger in (d). i | $m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$ | $m_{5}$ | $m_{6}$
---|---|---|---|---|---|---
$p_{i}$ | 0.988 | 0.980 | 0.975 | 0.979 | 0.971 | 0
$q_{i}$ | 0.988 | 0.993 | 0.994 | 0.998 | 0.996 | 0.89
Table 2:  Sensitivity and specificity parameters $p_{i}$ and $q_{i}$ for the
segmentations $m_{1}$ to $m_{6}$ displayed in Fig.8.

In Fig.9, we show the evolution of the active contour from the initial
contours (bubbles) given in Fig.9.a. One intermediate contour is given
Fig.9.b, and the final mutual shape is shown in Fig.9.c. The mutual shape
provided in Fig.9.d provides an interesting result for segmentation fusion
that takes benefit of the different segmentation entries while being robust to
the outlier shape. The evolution of the active contour displayed in Fig.9
shows that the initial shape evolves correctly towards the boundaries of the
object of interest. The sensitivity and specificity parameters are computed
together with this reference shape and provided in Table 2. These parameters
provide an interesting classification of the different segmentation without
any given reference. The mask $m_{2}$ seems to be the best segmentation
regarding with the estimated consensus and $m_{6}$ clearly appears as an
outlier. However, such a classification is dependent on the choice of the
different segmentation entries. We can however conclude that the mutual shape
is robust to the introduction of an outlier segmentation in the entry sequence
and that the outlier is clearly detected at the end of the process through the
low values of its corresponding $p_{i}$ ($p_{i}=0$ in Table 2).

Concerning the computational cost, for this example, the estimation of the
mutual shape takes $37s$ using Intel-based CPU @ 2.70GHz. The size of the
image is $321*481$.

###  7.2 Application to text segmentation in old parchments

A second real application of our algorithm is dedicated to the fusion of
different segmentations of the text in old parchments.

As a first example, we propose to combine different basic binarization methods
using the mutual shape in order to construct a better segmentation. Let us
consider for example the original image given in Fig.10 where the object of
interest is the whole text. The input masks are obtained using classical
binarization techniques provided by the library of image processing Pandore
[31]. The techniques used are namely "pmassbinarization” (based on a
percentage of pixels, mask 1 and 2), “pcorrelationbinarization” (maximization
of the correlation between two classes, mask 3), “pvariancebinarization”
(maximization of the interclass and intraclass distance, mask 4),
“pniblackbinarization” (based on an adaptive binarization technique described
in [36] mask 5) and “padaptativemeanbinarization” (based on the analysis of
the mean value of the intensities on a sliding window mask 6) . The
corresponding masks (shown in Fig.11) are used as segmentation inputs of our
mutual shape algorithm.

Figure 10: An original image from an old manuscript from Gallica (Gallica is the online numerical library of the BNF (National French Library)). |  |
---|---|---
(a) mask 1 | (b) mask 2 | (c) mask 3
|  |
(d) mask 4 | (e) mask 5 | (f) mask 6
Figure 11:  The different segmentation masks of the text (a,b,c,d,e,f).

The mutual shape is then computed using active contours (Fig.12). The initial
contour is chosen as a set of little circles currently named as “bubbles” in
the framework of active contours. The text is well segmented as displayed in
Fig.12.c showing the potential application of this method to build a consensus
segmentation from a set of different simple binarization techniques not
necesseraly all well chosen and composed of a set of pixels that is not
connected. This example also shows that our algorithm is able to handle a
shape composed of different separated components.

|
---|---
(a) Initial contour | (b) Iteration 50

---
(c) Final contour

---
(d) Final segmentation

Figure 12: Evolution of the active contour (in black) using the mutual shape
(evolution equation (25) with $\lambda=10$). The active contour is in black.
In the first image (a), the initial contour corresponds to a set of bubbles.
An Intermediate result obtained from $50$ iterations is displayed in image (b)
and the final contour in (c). The segmented object is displayed in (d) (the
background is in black color).

Let us now take another example of applicability of the mutual shape for
segmentation fusion and evaluation of different methods of text binarization
taken from the DIBCO database 2013 [33]. Let us consider for example the
original image given in Fig.13.a where the object of interest is the whole
text and for which we have the reference segmentation given in Fig.13.b.

|
---|---
(a) | (b)
Figure 13: An original image from an old manuscript from the DIBCO database
(a) and its corresponding reference binarization.

The input masks corresponds to different algorithms tested during this
challenge and are all available in the database. They are given in Fig.14.

|
---|---
(a) mask 1 | (b) mask 2
|
(c) mask 3 | (c) mask 4
|
(e) mask 5 | (f) mask 6
|
(g) mask 7 |
Figure 14:  The different segmentation masks of the text (a,b,c,d,e,f,g) taken
from the 2013 DIBCO database.

For all these masks, we can compute the Dice Coefficient with the reference
segmentation. The different values are given in Table.3.

i | $m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$ | $m_{5}$ | $m_{6}$ | $m_{7}$
---|---|---|---|---|---|---|---
$DC$ | 0.54 | 0.89 | 0.95 | 0.87 | 0.91 | 0.76 | 0.91
Table 3: Computation of the Dice Coefficient with the reference segmentation
(Fig.13.b) for the segmentations $m_{1}$ to $m_{7}$ (Fig.14).

Let us now compute the mutual shape and compare the quality of the obtained
result to the reference segmentation. The obtained mutual shape (final contour
and the associated mask) are given in Fig.15. For this mutual shape, we find
$DC=0.93$ which outperforms the DC coefficient of all the different masks in
entry except the mask $m_{3}$. The resulting mutual shape is then interesting
for an intelligent fusion of different segmentation results. Let us now
compare the ranking obtained using the mutual shape algorithm through the
joint computation of the $p_{i}$ and $q_{i}$ coefficients. The different
values of $p_{i}$ and $q_{i}$ are given in Table.4 and allow us to rank the
different segmentation methods as follows (from the best one to the worst one
according to the sum of $p_{i}$ and $q_{i}$) :
$m_{3}$,$m_{5}$,$m_{2}$,$m_{7}$,$m_{6}$,$m_{4}$,$m_{1}$. The ranking obtained
using the reference mask and the DC coefficients leads to
:$m_{3}$,$m_{7}$,$m_{5}$,$m_{2}$,$m_{4}$,$m_{6}$,$m_{1}$. We can observe that
the ranking is the same for the first and the last mask. There are some
difference of ranking between comparable masks such as $m_{7}$, $m_{2}$ and
$m_{5}$. The mask $m_{4}$ corresponds to an under-segmentation and the mask
$m_{6}$ to an over-segmentation which explains their places in the end of the
ranking.

---
(a) Final contour
(b) Resulting mask
Figure 15: The mutual shape (final contour and resulting mask) computed using the masks $m_{1}$ to $m_{7}$ of the DIBCO image PR08. i | $m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$ | $m_{5}$ | $m_{6}$ | $m_{7}$
---|---|---|---|---|---|---|---
$p_{i}$ | 0.973 | 0.975 | 0.984 | 0.872 | 0.982 | 0.992 | 0.969
$q_{i}$ | 0.862 | 0.985 | 0.993 | 0.997 | 0.989 | 0.949 | 0.990
Table 4:  Sensitivity and specificity parameters $p_{i}$ and $q_{i}$ for the
segmentations $m_{1}$ to $m_{7}$ displayed in Fig.14.

In this last example, we show that the mutual shape leads to an interesting
segmentation result by performing an intelligent fusion of different
segmentation entries. The obtained ranking is interesting but can be different
to the ranking performed using a reference mask and the DC coefficient.

###  7.3 Application to segmentation of cardiac magnetic resonance images

The segmentation of cardiac structures is an active research field in all
medical modalities [16], where expert performance is still higher than image
segmentation algorithms performance. As experts segmentations can vary, it was
proposed to use STAPLE algorithm to define a consensus segmentation between
different experts [39]. Furthermore, to reduce the drawbacks of each specific
image segmentation algorithm, it was proposed to take advantage of the results
of different segmentation algorithms and the interest of combining different
segmentation results using STAPLE was shown [26]. When compared to individual
methods, using these combined segmentations provided better estimates of the
clinical parameters of interest; this was demonstrated by a supervised
approach using experts delineations and a non supervised evaluation approach
described in [25]. In this specific context, the estimation of a mutual shape
was tested for the non supervised evaluation and the fusion of different
segmentation methods of the left ventricular cavity from cardiac cine magnetic
resonance images [20, 19]. For instance, the excellent behavior of mutual
shape towards outliers was demonstrated. In this section, we propose a first
comparison between the mutual shape approach and STAPLE. At the difference of
mutual shape, the STAPLE algorithm does not introduce any regularization term
and can thus provide unsmoothed results, which are not relevant on a
physiological basis.

The segmentation entries are selected inside a database that contains the
results obtained by three experts and different algorithms [9, 11, 24, 37,
10]. The corresponding contours are displayed in Fig.16 and Fig.17. In this
specific example, the endocardium is not well delimited by the automated
algorithms, due to the presence of the aortic root, which leads to very
different segmentations.

|  |
---|---|---
$Exp_{1}$ | $Exp_{2}$ | $Exp_{3}$
Figure 16: Segmentation methods of the left ventricle provided by three experts $Exp_{1}$ to $Exp_{3}$.  |  |  |  |
---|---|---|---|---
$m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$ | $m_{5}$
Figure 17: Segmentation methods of the left ventricle provided by $5$ automated algorithms $m_{1}$ to $m_{5}$.  |  |  |  |
---|---|---|---|---
(a) | (b) | (c) | (d) | (e)
Figure 18: Estimation of different consensus estimates using the contours given by the three experts $Exp_{1}$ to $Exp_{3}$ (Fig.16) using STAPLE algorithm (a), SD approach (b), and mutual shape (c). Filled masks correspond to STAPLE (d) and mutual shape (e) |  |  |
---|---|---|---
(a) | (b) | (c) | (d)
Figure 19: Estimation of different consensus estimates using the masks $m_{1}$ to $m_{5}$ (Fig.17) using STAPLE algorithm (a), SD approach (b), and mutual shape (c), contour delineated by $Exp_{2}$ (d). $i$ | $Exp_{1}$ | $Exp_{2}$ | $Exp_{3}$
---|---|---|---
$p_{i}$ | 0.999 | 0.999 | 0\. 999
$q_{i}$ | 0.889 | 0.973 | 0.853
Table 5: Sensitivity and specificity parameters $p_{i}$ and $q_{i}$ for the segmentations $Exp_{1}$ to $Exp_{3}$ (Fig.16). $i$ | $m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$ | $m_{5}$
---|---|---|---|---|---
$p_{i}$ | 0\. 884 | 0.715 | 0.956 | 0.714 | 0.787
$q_{i}$ | 0.847 | 0.999 | 0.784 | 0.999 | 0.999
Table 6: Sensitivity and specificity parameters $p_{i}$ and $q_{i}$ for the
segmentations $m_{1}$ to $m_{5}$ displayed in Fig.17.

Fig.18 shows the consensus segmentations estimated by STAPLE (a), SD (b), and
mutual shape (c), using the three experts entries. The regularization
parameter $\lambda$ was set equal to $100$. Table 5 shows that $Exp_{2}$
provided for this specific case the best contour, but this result (superiority
of $Exp_{2}$ ) was already reported elsewhere [26]. Filled masks of STAPLE
(Fig.18.d) and mutual shape (Fig.18.e) clearly demonstrate that STAPLE does
not necessarily provide smoothed contours at the difference of mutual shape.

Finally, Fig.19 shows the consensus segmentations, estimated by STAPLE (a), SD
(b), and mutual shape (c), using the five segmentation entries of the
automated algorithms ($m_{1}$ to $m_{5}$). The regularization parameter
$\lambda$ was set equal to $10$. The parameters $p_{i}$ and $q_{i}$ are
estimated jointly with the mutual shape (see Table 6). At the difference of
STAPLE, the contour provided by mutual shape is smooth. However both methods
provide results that are quite different from the expert entries. This
difficult case shows that mutual shape and STAPLE both depend on the accuracy
of the segmentation entries.

##  8 Conclusion

In this work, we search for a mutual shape that minimizes the sum of joint
entropies while maximizing the sum of mutual information between each entry
shape and the unknown reference shape. We give a geometrical interpretation of
this criterion using area measures. The optimization is performed using active
contours by computing a shape gradient and the associated evolution equation.
Shape derivatives are computed and detailed for the given criterion. Our
theoretical formalism is valid for 2D slices or 3D images. The main
contribution of this paper lies in the proposition of a theoretical criterion
for the estimation of a consensus shape both for segmentation fusion and
evaluation without reference.

Some experimental results are provided on both synthetic images and real
images for segmentation fusion and evaluation. Indeed, the proposed mutual
shape is able to build a consensus shape from a set of different segmentations
of the same object and can then be used as an intelligent fusion of different
segmentation entries. Moreover, the algorithm estimates jointly the
sensibility and specificity parameters (as first proposed by Warfield et al
[41]). It then provides an evaluation or ranking of the proposed segmentation
methods on the basis of these two parameters.

The experimental results performed on a synthetic image allow to better
understand the difference between the mutual shape, a simple union, an average
shape and a simple majority voting method. Moreover, the estimated mutual
shape is robust to very aberrant outliers thanks to the joint minimization of
the joint entropy with the maximization of the mutual information.

We then propose some tests on real images for different applications. A first
test is performed on a real color image and we show that the estimated mutual
shape is robust and accurate. The second type of application concerns the
segmentation of the text in old manuscripts from a set of simple segmentation
entries. The results obtained for this application bring new opportunitues in
the field of segmentation by demonstrating that an improved segmentation
method may be designed by taking benefit of the fusion of several simple
segmentation methods. The intelligent fusion of simple algorithms can probably
lead to a new powerful segmentation algorithm. For this part, the choice of
the segmentation entries is still an open issue but one may think of a
learning phase for an interactive choice of the different segmentation methods
chosen to build the consensus. We may call this new kind of process :
segmentation by intelligent voting using a consensus shape. The last
application, devoted to medical images, also shows that the mutual shape can
be useful for the evaluation of different segmentation methods without any
reference. However, such a classification may be a little different than the
one obtained using the true reference shape. Indeed it corresponds to a
classification regarding with the most preponderant shapes in the set of
entries and not to an absolute classification with a fixed reference. The
classification is clearly dependent on the choice of the different
segmentation methods in entry. This unsupervised evaluation process may
however be useful to detect abnormal segmentation methods in a set of
different segmentation entries.

One perspective of this work may also concern the addition of prior terms
inside the variational criterion. For example, a shape prior can be
interesting when the objective is to segment the left ventricular cavity. If
the objective is different, some other prior shapes may be added (such as the
homogeneity of the inside region for example, the gradient, or the target
color for color segmentation). Our mathematical framework seems well adapted
for this purpose since other information may be easily and rigorously
introduced in the criterion to minimize.

## Acknowledgements

The medical application took place in a larger project named MediEval
supported by the GdR 2647 Stic Santé (CNRS-INSERM). We thank all the partners
of this project for providing us the different segmentations for the
evaluation part.

## Appendix

### Proof of Theorem 3

First of all, we compute the shape derivatives of the probabilities $p_{i}$
and $q_{i}$ which depend on the domain. By applying the theorem 1, we find :

$\displaystyle<p_{i}^{\prime}(\mu),{\mathbf{V}}>=-\frac{1}{|\mu|}\int_{\partial\mu}K(d_{i}({\mathbf{x}})-1)({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$
$\displaystyle+\frac{\int_{\partial\mu}({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}}{|\mu|^{2}}\int_{\mu}K(d_{i}({\mathbf{x}})-1)d{\mathbf{x}}$

which reduces to :

$\displaystyle<p_{i}^{\prime}(\mu),{\mathbf{V}}>=\frac{1}{|\mu|}\int_{\partial\mu}(p_{i}-K(d_{i}({\mathbf{x}})-1))({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$
(29)

In the same way, we compute the shape derivative of $q_{i}$ :

$\displaystyle<q_{i}^{\prime}(\mu),{\mathbf{V}}>=\frac{1}{|\overline{\mu}|}\int_{\partial\overline{\mu}}(-q_{i}+K(d_{i}({\mathbf{x}})))({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$
(30)

We can also compute the shape derivatives of $\mu$ and $\overline{\mu}$:

$\displaystyle<|\mu|^{\prime},{\mathbf{V}}>=-\int_{\partial\mu}({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$
$\displaystyle<|\overline{\mu}|^{\prime},{\mathbf{V}}>=\int_{\partial\mu}({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$

Let us denote $\varphi(p)=p\log(p)+(1-p)\log(1-p)$, the conditional entropy
then becomes :

$H(D_{i}/T)=-\Big{[}\frac{|\mu|}{|\Omega|}\varphi(p_{i})+\frac{|\overline{\mu}|}{|\Omega|}\varphi(q_{i})\Big{]}$

By using chain derivation rules, we find :

$\displaystyle<H(D_{i}/T)^{\prime},{\mathbf{V}}>$ $\displaystyle=$
$\displaystyle-\frac{1}{|\Omega|}\Big{[}|\mu|<p_{i}^{\prime},{\mathbf{V}}>\varphi^{\prime}(p_{i})$
$\displaystyle+$
$\displaystyle|\overline{\mu}|<q_{i}^{\prime},{\mathbf{V}}>\varphi^{\prime}(q_{i})+<|\mu|^{\prime},{\mathbf{V}}>\varphi(p_{i})$
$\displaystyle+$
$\displaystyle<|\overline{\mu}|^{\prime},{\mathbf{V}}>\varphi(q_{i})\Big{]}$

where $\varphi^{\prime}(p)=\log(p)-\log(1-p)$. Replacing the shape derivatives
by the previous formulas, we find Theorem 3.

### Proof of Theorem 4

First of all, we compute the shape derivatives of the joint probabilities
$p(d_{i},t)$ which depend on the domain. Using theorem 1, we find for a=1 or
a=0 :

$<p(d_{i}=a,t=0)^{\prime},{\mathbf{V}}>=-\frac{1}{|\Omega|}\int_{\partial\mu}K(d_{i}({\mathbf{x}})-a)({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$

and

$<p(d_{i}=a,t=1)^{\prime},{\mathbf{V}}>=\frac{1}{|\Omega|}\int_{\partial\mu}K(d_{i}({\mathbf{x}})-a)({\mathbf{V}}\cdot{\mathbf{N}})d\mathbf{a}$

Let denote $\Psi(p)=p\log p$, the joint entropy then becomes :

$\displaystyle
H(D_{i},T)=-\Psi\left((p(d_{i}=0,t=0)\right)-\Psi\left(p(d_{i}=1,t=0)\right)$
$\displaystyle-\Psi\left(p(d_{i}=0,t=1)\right)-\Psi\left(p(d_{i}=1,t=1)\right)$

By using chain derivation rules, we find :

$\displaystyle<H(D_{i},T)^{\prime},{\mathbf{V}}>=-<p(d_{i}=1,t=1)^{\prime},{\mathbf{V}}>\Psi^{\prime}(p(d_{i}=1,t=1))$
$\displaystyle-<p(d_{i}=0,t=1)^{\prime},{\mathbf{V}}>\Psi^{\prime}(p(d_{i}=0,t=1))$
$\displaystyle-<p(d_{i}=1,t=0)^{\prime},{\mathbf{V}}>\Psi^{\prime}(p(d_{i}=1,t=0))$
$\displaystyle-<p(d_{i}=0,t=0)^{\prime},{\mathbf{V}}>\Psi^{\prime}(p(d_{i}=0,t=0))$

where $\Psi^{\prime}(p)=\log(p)+1$. Replacing the shape derivatives by the
previous formulas, we find Theorem 4.

## References

  * [1] G. Aubert, M. Barlaud, O. Faugeras, and S. Jehan-Besson.  Image segmentation using active contours : Calculus of variations or shape gradients.  SIAM Journal on Applied Mathematics, 63(6):2128–2154, 2003.
  * [2] B. Berkels, G. Linkmann, and M. Rumpf.  An SL(2) invariant shape median.  Journal of Mathematical Imaging and Vision, 37(2):85–97, 2010.
  * [3] X. Bresson, S. Esedoglu, P. Vandergheynst, J. Thiran, and S. Osher.  Fast global minimization of the active contour/snake model.  Journal of Mathematical Imaging and Vision, 28(2):151–167, 2007\.
  * [4] T. F. Chan and L. A. Vese.  Active contour without edges.  IEEE Transactions on Image Processing, 10(2):266–277, 2001.
  * [5] T.F. Chan, B.Y. Sandberg, and L.A. Vese.  Active contours without edges for vector-valued images.  Journal of Visual Communication and Image Representation, 11(2):130–141, 2000.
  * [6] G. Charpiat, O. Faugeras, and R. Keriven.  Shape metrics, warping and statistics.  In International Conference on Image Processing, volume 2, pages 627–630, 2003.
  * [7] G. Charpiat, O. Faugeras, and R. Keriven.  Approximations of shape metrics and application to shape warping and empirical shape statistics.  Foundations of Computational Mathematics, 5(1):1–58, 2005.
  * [8] O. Commowick, A. Akhondi-Asl, and S.K. Warfield.  Estimating a reference standard segmentation with spatially varying performance parameters.  IEEE Transactions on Medical Imaging, 31(8):1593–1606, 2012.
  * [9] C. Constantinides, R. El Berbari, A. de Cesare, Y. Chenoune, E. Roullot, A. Herment, E. Mousseaux, and F. Frouin.  Development and evaluation of an algorithm for the automated segmentation of the left and right ventricles on cine MRI.  IRBM, 30(4):188–191, 2009.
  * [10] C. Constantinides, E. Roullot, M. Lefort, and F. Frouin.  Fully automated segmentation of the left ventricle applied to cine MR images: description and results on a database of 45 subjects.  In IEEE Engineering in Medicine and Biology Society, pages 3207–3210, 2012.
  * [11] J. Cousty, L. Najman, M. Couprie, S. Clément-Guinaudeau, T. Goissen, and J. Garot.  Segmentation of 4D cardiac MRI: automated method based on spatio-temporal watershed cuts.  Image and Vision Computing, 28(8):1229–1243, 2010.
  * [12] T.M. Cover and J.A. Thomas.  Elements of Information Theory.  Wiley-Interscience, 1991.
  * [13] D. Cremers, T. Kohlberger, and C. Schnörr.  Shape statistics in kernel space for variational image segmentation.  Pattern Recognition, 36(2):1929–1943, 2003.
  * [14] M.C. Delfour and J.P. Zolésio.  Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization.  Advances in design and control. Society for Industrial and Applied Mathematics (SIAM), 2001.
  * [15] R. Duda and P. Hart.  Pattern Classification and Scene Analysis.  John Wiley & Sons, Inc., 1973.
  * [16] J. Fleureau, M. Garreau, D. Boulmier, C. Leclercq, and A. Hernandez.  3D multi-object segmentation of cardiac MSCT imaging by using a multi-agent approach.  IRBM, 30(3):104–113, 2009.
  * [17] A. Foulonneau, P. Charbonnier, and F. Heitz.  Affine-invariant geometric shape priors for region-based active contours.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(8):1352–1357, 2006.
  * [18] A. Herbulot, S. Jehan-Besson, S. Duffner, M. Barlaud, and G. Aubert.  Segmentation of vectorial image features using shape gradients and information measures.  Journal of Mathematical Imaging and Vision, 25(3):365–386, 2006\.
  * [19] S. Jehan-Besson, C. Tilmant, A. de Cesare, A. Lalande, A. Cochet, J. Cousty, J. Lebenberg, M. Lefort, P. Clarysse, R. Clouard, L. Najman, L. Sarry, F. Frouin, and M. Garreau.  A mutual reference shape based on information theory.  In IEEE International Conference on Image Processing, pages 887–891, 2014.
  * [20] S. Jehan-Besson, C. Tilmant, A. De de Cesare, F. Frouin, L. Najman, A. Lalande, L. Sarry, C. Casta, P. Clarysse, C. Constantidinès, J. Cousty, M. Lefort, A. Cochet, and M. Garreau.  Estimation d’une forme mutuelle pour l’évaluation de la segmentation en imagerie cardiaque.  In GRETSI - Traitement du Signal et des Images, 2011.
  * [21] M. Kass, A. Witkin, and D. Terzopoulos.  Snakes : Active contour models.  International Journal of Computer Vision, 1:321–332, 1988.
  * [22] D.G. Kendall.  Shape manifolds, procrustean metrics, and complex projective spaces.  Bulletin of the London Mathematical Society, 16(2):81–121, 1984\.
  * [23] J. Kim, J.W. Fisher, A. Yezzi, M. Cetin, and A.S. Willsky.  Nonparametric methods for image segmentation using information theory and curve evolution.  In IEEE International Conference on Image Processing, pages 797–800, 2002.
  * [24] A. Lalande, N. Salve, A. Comte, M.-C. Jaulent, L. Legrand, P.M. Walker, Y. Cottin, J.E. Wolf, and F. Brunotte.  Left ventricular ejection fraction calculation from automatically selected and processed diastolic and systolic frames in short axis cine-MRI.  Journal of Cardiovascular Magnetic Resonance, 6(4):817–827, 2004\.
  * [25] J. Lebenberg, I. Buvat, A. Lalande, P. Clarysse, C. Casta, A. Cochet, C. Constantinidès, J. Cousty, A. de Cesare, S. Jehan-Besson, M. Lefort, L. Najman, E. Roullot, L. Sarry, C. Tilmant, M. Garreau, and F. Frouin.  Non supervised ranking of different segmentation approaches: application to the estimation of the left ventricular ejection fraction from cardiac cine MRI sequences.  IEEE Transactions on Medical Imaging, 31(8):1651–1660, 2012.
  * [26] J. Lebenberg, A. Lalande, P. Clarysse, I. Buvat, C. Casta, A. Cochet, C. Constantinidès, J. Cousty, A. de Cesare, S. Jehan-Besson, M. Lefort, L. Najman, E. Roullot, L. Sarry, C. Tilmant, F. Frouin, and M. Garreau.  Improved estimation of cardiac function parameters using a combination of independent automated segmentation results in cardiovascular magnetic resonance imaging.  PLOS ONE, 10(8):e0135715, 2015.
  * [27] F. Lecellier, S. Jehan-Besson, J. Fadili, G. Aubert, M. Revenu, and E. Saloux.  Region-based active contours with noise and shape priors.  In IEEE International Conference on Image Processing, volume 1, pages 1649–1652, 2006.
  * [28] M. Leventon, E. Grimson, and O. Faugeras.  Statistical shape influence in geodesic active contours.  In IEEE International Conference on Image Processing, pages 316–323, 2000.
  * [29] D. Martin, C. Fowlkes, D. Tal, and J. Malik.  A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics.  In Proceedings Eighth International Conference on Computer Vision, volume 2, pages 416–423, 2001.
  * [30] S. Osher and J.A. Sethian.  Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations.  Journal of Computational Physics, 79(1):12–49, 1988.
  * [31] Pandore: A library of image processing operators (Version 6.6). [Software]. Greyc Laboratory.  https://clouard.users.greyc.fr/Pandore, 2013.
  * [32] N. Paragios.  A level set approach for shape-driven segmentation and tracking of the left ventricle.  IEEE Transactions on Medical Imaging, 22(6):773–776, 2003.
  * [33] I. Pratikakis, B. Gatos, and K. Ntirogiannis.  ICDAR 2013 document image binarization contest (DIBCO 2013).  In Proceedings of the 2013 12th International Conference on Document Analysis and Recognition, pages 1471–1476, Washington, DC, USA, 2013\. IEEE Computer Society.
  * [34] F. Precioso, M. Barlaud, T. Blu, and M. Unser.  Robust real-time segmentation of images and videos using a smooth-spline snake-based algorithm.  IEEE Transactions on Image Processing, 14(7):910–924, 2005.
  * [35] F. M. Reza.  An Introduction to Information Theory, pages 106–108.  McGraw-Hill, 1994.
  * [36] J. Sauvola and M. Pietikainen.  Adaptive document image binarization.  Pattern Recognition, 33(2):225–236, 2000.
  * [37] J. Schaerer, C. Casta, J. Pousin, and P. Clarysse.  A dynamic elastic model for segmentation and tracking of the heart in MR image sequences.  Medical Image Analysis, 14(6):738–749, 2010.
  * [38] J. Sokolowski and J.P. Zolésio.  Introduction to shape optimization. Shape sensitivity analysis. Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1992.
  * [39] A. Suinesiaputra, B.R. Cowan, A.O. Al-Agamy, M.A. Elattar, N. Ayache, A.S. Fahmy, A.M. Khalifa, P. Medrano-Gracia, M.P. Jolly, A.H. Kadish, D.C. Lee, J. Margeta, S.K. Warfield, and A.A. Young.  A collaborative resource to build consensus for automated left ventricular segmentation of cardiac MR images.  Medical Image Analysis, 18(1):50–62, 2014.
  * [40] S. Velasco-Forero and J. Angulo.  Statistical shape modeling using morphological representations.  In IEEE International Conference on Image Processing, pages 3537–3540, 2010.
  * [41] S.K. Warfield, K. H. Zou, and W. M. Wells III.  Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation.  IEEE Transactions on Medical Imaging, 23(7):903–921, 2004.
  * [42] R. W. Yeung.  A new outlook on shannon’s information measures.  IEEE Transactions on Information Theory, 37(3):466–474, 1991.
  * [43] A.J. Yezzi and S. Soatto.  Deformotion: Deforming motion, shape average and the joint registration and approximation of structures in images.  International Journal of Computer Vision, 53(2):153–167, 2003.


# Estimation of Radio Refractivity from a Decade Satellite Derived
Meteorological Data for West Africa

O. J. Abimbola<EMAIL_ADDRESS>E. Bada  A. O.
Falaiye  Y. M. Sukam  S. Muhammad  Department of Physics, Federal University
of Lafia, Nigeria  Department of Physics, University of Ilorin, Nigeria

###### Abstract

Radio refractivity, which is the bending of radio signal as it propagates
through media, is very important in works involving terrestrial atmospheric
electromagnetic propagation such as point-to-point microwave communication,
terrestrial radio and television broadcast, mobile communication system, and
so on. This study has focused on the West African region where it was found
that the region has a refractivity which varies exponentially with height and
from the coast towards the desert. Generally, refractivity gradient was found
to range between -46.48 and -29.51 N-units/km (k-factor value of between 1.23
and 1.42) across the region, splitting the between sub- and super-refraction.
The variation in refraction type was found to follow seasonal pattern across
the West African region, with sub-refraction dominating during dry season and
super-refraction dominating most part in the coastal area during the wet
season.

###### keywords:

radio refractivity, West Africa, refractivity gradient, k-factor

††journal: arXiv

##  1 Introduction

The behavior of radio waves in the tropospheric layer of the Earth’s
atmosphere is very important in this modern age that is highly influenced by
radio communications ranging from mobile telephoning through terrestrial
digital broadcasting to the propagation of satellite radio signal through the
troposphere.

Propagation of radio waves in the inhomogeneous troposphere leads to the
refraction of the radio waves with subsequent refractive index n, given by

$n=\frac{speed~{}of~{}EM-waves~{}in~{}air}{speed~{}of~{}EM-
waves~{}in~{}vacuum}$ (1)

The refractive index n, usually have a very small value for the troposphere
(1.0000 to 1.0000004445). In order to have an easy to work with value, n is
usually transformed to the refractivity N by [1]

$N=\left(n-1\right)\times 10^{6}$ (2)

Since n is a unit-less quantity, N is also expected to be unit-less,
nevertheless, N is usually given in the unit of “N-unit”. In the troposphere,
N ranges between 200 to 400 N-units.

According to the International Telecommunication Union recommendation [2], the
radio refractivity index N, can be given as N(p, e, T), where p is the total
atmospheric pressure, e is the vapour pressure and T is the ambient
temperature:

$N=77.6\frac{p}{T}+3.75\times 10^{5}\frac{e}{T^{2}}$ (3)

where, p and e are measured in hPa and T is measured in Kelvin. As is obvious
from Equation 3, atmospheric radio refractivity depends significantly on the
atmospheric moisture content. It should be noted that Equation 3 is valid for
radio frequencies up to $100~{}GHz$ [2, 3, 4].

In grouping the radio refractivity into groups such as normal/standard
refraction, sub-refraction, super-refraction and ducting, a quantity known as
the effective earth radius factor or k-factor, is employed. The k-factor is
given as

$k=\frac{1}{1+\left(dN/dh\right)/157}$ (4)

here, $dN\textfractionsolidus dh=G$ which is known as the refractivity
gradient [4, 5, 6]. G is a quantity which is very useful in describing radio
refractivity at the first $1~{}km$ of the troposphere, from the earth’s
surface [2]: if $G>-40$ then refraction is sub-refraction; if $G<-40$ then
refraction is super-refraction while if $G<-157$ the refraction is termed
ducting [7]. In terms of earth radius factor, standard refractivity occurs at
$k=1.33$; sub-refraction, that is, bending of the radio waves away from the
earth’s surface, occurs at $1.33<k<\infty$; bending of the radio waves towards
the earth, known as super-refraction, occurs at $0<k<1.33$ while, ducting
occurs at $-\infty<k<0$ [4, 8].

Several investigators have worked on the study of radio refractivity using
localized data from radiosonde and masts. Davis et al. [9] used ground-based
microwave radiometer to study the local spatio-temporal variation of the radio
refractivity in west coast of Sweden, they found an average surface wet-
refractivity horizontal gradients of between 0.1 – 1.0 $N-units/km$ by
assuming exponential profile for the wet-refractivity horizontal gradient,
however they noted that their model will be inadequate for period of time
longer than 30 $min$. Focusing on some locations within the North Central
Nigeria, Emetere et al. [10] studied the effects of aerosol on the radio
refractivity. Ojo et al. [11], studied geoclimatic variation of radio
refractivity across Nigeria where map of k-factor for Nigeria was presented
and they generally found that k-factor always exceed the value of 1.33 in
Nigeria.

It is our aim, in this study, to use a ten-year satellite derived data to
describe both horizontal and vertical distribution of radio refractivity in
West Africa together with the description of refractivity gradient and the
earth radius factor.

##  2 Methodology

West African region under consideration is between $18~{}^{o}W$ and
$18~{}^{o}E$ longitudes. This region is a low latitude region between latitude
$3~{}^{o}N$ and $18~{}^{o}N$. The study area generally has two distinct
seasons: the wet season and the dry season. The wet season is a period of high
atmospheric humidity and rainfall, the period vary across the region, at the
coast the wet season could occur within the months of April to October while
at places closer to the Sahara Desert the wet season could be from June to
September. The remaining part of the year is marked by low atmospheric
humidity and practically no rainfall, the early part of the dry season, within
November, December and January, there is heavy dust lifting from the Sahara
Desert and blown across the region, this dust movement period is called
Harmattan period: this seasonal variation of atmospheric moisture
significantly affects radio refractivity.

Figure 1: Map of West Africa showing the stations used for this study, marked
in red dots

Ambient temperature and specific humidity data were obtained from the archive
of the Climate Monitoring Satellite Application Facility (CM-SAF) of the
European Center for Medium Range Weather Forecast (ECMWF) at
https://wui.cmsaf.eu/safira/action while the latitudinal radio refractivity
was obtained from the Radio Occultation Meteorology Satellite Application
Facility (ROM SAF) at https://romsaf.org. CM-SAF data was obtained at six
different vertical pressure levels: 200 hPa, 300 hPa, 500 hPa, 700 hPa, 850
hPa and 1000 hPa, as a daily average for the years 2007 to 2016. The data from
CM-SAF was obtained in netCDF format and python script together with Microsoft
Excel software was used to extract the data while Equation 5 was used to
calculate the water vapour partial pressure e

$e=\frac{0.001pq}{0.622}$ (5)

where, p = total atmospheric pressure; q = specific humidity. All geo-plots
were done using QGIS freeware.

##  3 Results and Discussion

Figure 2 shows the spatial distribution of the average radio refractivity for
the years 2006 to 2016 at four different pressure levels. Figure 3 shows the
latitudinal vertical variation of radio refractivity during the year 2018 for
the months of January, April, August and December: generally in West Africa,
December and January are within the dry very low humidity period of the year
while April is the beginning of the wet humid season and August is about the
peak of the wet humid season. Exponential vertical variation of radio
refractivity at latitudes $2.5~{}^{o}N,7.5~{}^{o}N,12.5~{}^{o}N,17.5~{}^{o}N$
and $22.5~{}^{o}N$, is clearly visible on Figure 4.

Figure 2: Spatial distribution of radio refractivity at different pressure
level within the West African region: (A) at 850 hPa pressure level (that is,
around 1.54 km amsl), (B) at 700 hPa pressure level (that is, around 3.18 km
amsl), (C) at 500 hPa pressure level (that is, around 9.89 km amsl), (D) at
200 hPa pressure level (that is, around 12.14 km amsl).

At 850 hPa, that is about 1.54 km above mean sea level, as presented on Figure
2A, the radio refractivity was found to range between 306.67 and 257.69
N-units with a range of 48.98 N-units, the radio refractivity could be
observed to be higher at the coast and gradually decreases as we move to a
dryer locations closer to the Sahara Desert, this is as a result of the
dependence of radio refractivity on the atmospheric moisture burden. This
variation in radio refractivity is replicated across other pressure levels but
at 200 hPa, that is, about 12.14 km above mean sea level, where the radio
refractivity ranges between 71.37 and 71.14 N-units with a range of 0.23
N-units, the radio refractivity was found to be slightly higher east of the
prime meridian than to the west of the prime meridian: the coastal bound at
the western part is more than at the eastern part of the West African region
resulting in slightly more moisture at this altitude on the western part than
at the eastern part.

Figure 3: Vertical distribution of mean monthly radio refractivity, along the
latitudes for the months of January, April, August and December in the year
2018.

The vertical variation of the radio refractivity along the latitudes that
includes the West Africa region (i.e., 2.5, 7.5, 12.5, 17.5 and 22.5
${}^{o}N$) is as shown on Figure 3 and the mathematical models for the
vertical variations is as shown on Table 1. From Figure 4, it could be
observed that the vertical variations are exponential in nature which is what
is obtained and shown on Table 1. The surface value of radio refractivity was
found to range between 403 and 420 N-units while the radio refractivity scale
height was found to be around 6.6 km. It should be noted that the models shown
on Table 1 are for the latitudes and not regional, besides the table uses ROM-
SAF data.

Figure 4: Vertical variation of the radio refractivity at several latitudes (ROM-SAF data). Table 1: Mathematical models describing the vertical variation of the radio refractivity at different latitudes (ROM-SAF data). Latitude (${}^{o}N$) | Model | Coefficient of
---|---|---
|  | Determination ($R^{2}$)
2.5 | $N=420e^{-\left(H/6.6\right)}$ | 0.998
7.5 | $N=420e^{-\left(H/6.6\right)}$ | 0.998
12.5 | $N=415e^{-\left(H/6.6\right)}$ | 0.998
17.5 | $N=415e^{-\left(H/6.6\right)}$ | 0.997
22.5 | $N=403e^{-\left(H/6.6\right)}$ | 0.998
Table 2: Mathematical models describing the vertical variation of the radio refractivity at different locations within the West African region (CM-SAF data). Latitude | Longitude | Model | $R^{2}$
---|---|---|---
16$~{}^{o}N$ | 10$~{}^{o}E$ | $N=313e^{-\left(H/8.33\right)}$ | 0.999
6$~{}^{o}N$ | 5$~{}^{o}E$ | $N=363e^{-\left(H/7.69\right)}$ | 0.999
6$~{}^{o}N$ | 5$~{}^{o}W$ | $N=367e^{-\left(H/7.69\right)}$ | 0.998
16$~{}^{o}N$ | 10$~{}^{o}W$ | $N=326e^{-\left(H/8.33\right)}$ | 0.998

Using CM-SAF data, location based vertical models for the radio refractivity
were obtained for locations West and East of the prime meridian within the
West African region, as presented in Table 2. Within the considered locations,
the radio refractivity scale height was found to range between 7.28 and 8.19
km, while the surface refractivity ranges between 315 and 368 N-units, it is
obvious from both Table 1 and 2 that the vertical variations is exponential in
nature.

###  3.1 Refractivity Gradient and k-factor

A ten year average of the refractivity gradient across West Africa for the
first 1.5 km of the atmosphere, from the surface, is shown in Figure 5. From
Figure 5, it could be seen that two major types of refraction dominates the
atmosphere of West Africa, namely sub-refraction
($G>-40~{}N-units/km~{}or~{}\infty<k<1.33$) and super-refraction
($G<-40~{}N-units/km~{}or~{}0>k>1.33$). Figure 5B shows finer details than
Figure 5A, it could be observed that Figure 5B shows more area with super-
refraction.

Taking some countries as point of focus, for instance Nigeria, we could see,
on both Figure5A and B that the country is roughly divided into two parts, the
northern half is dominated by sub-refraction while the southern part is
dominated by super-refraction; this is in agreement with the work of Aremu et
al. [8], who worked in Ibadan, a southern Nigerian city and found G = -50.42
N-units/km as well as k = 1.53, also Adediji et al. [4] working in Akure, the
same southern part of Nigeria, found a two year average values of G and k to
be -52.8 N-units/km and 1.51 respectively. Fashade et al. [12], using a five
year satellite data found that Sokoto and Maiduguri, two cities in northern
Nigeria has average G values of -26.4 and -23.2 N-units/km respectively while
Abuja, Enugu, Lagos and Port Harcourt, cities in the middle belt (Abuja) and
southern Nigeria were found to have average G values of -40.2, -40.3, -47.0
and -43.0 N-units/km respectively. This result of Fashade et al. [12] agrees
very well with the result presented on Figure 5.

Figure 5: (A) West African refractivity gradient G map and (B) earth radius
factor k map, for the first 1.5 km of the atmosphere. Figure 6: Wet Season map
of (A) West African refractivity gradient G map and (B) earth radius factor k
map, for the first 1.5 km of the atmosphere.

If Figure 5B is used for study, Senegal, Gambia, Guinea-Bissau and the western
part of Guinea could be observed to be generally dominated by super-
refraction, this is also true of Benin, Togo, Ghana and eastern part of Cote
de Voire. The remaining countries of West Africa are then found to be
dominated by sub-refraction, the region with sub-refraction have values of
k-factor ranging between 1.23 and 1.32 (G ranges between -39.88 and -29.51
N-units/km) while the region of super-refraction have k-factor ranging between
1.34 and 1.42 (G ranges between -46.48 and -40.82 N-units/km). The regions of
sub- and super-refractions could be attributed to the atmospheric moisture
distribution in the West Africa region; the region of super-refraction could
be observed to be close to the Atlantic Ocean and hence have higher
atmospheric moisture burden while the region dominated by sub-refraction is a
region much closer to the Sahara Desert, hence with drier atmosphere.

###  3.2 Seasonal Variations of Refractivity Gradient and k-factor

Figure 6 shows the distribution of refractivity gradient and k-factor across
West Africa during the wet season months of May, June, July, August, September
and October. As could be observed on Figure 6B, the whole of Nigeria, Benin,
Togo and Ghana as well as major part of Cote de Voire, Senegal, Gambia, and
Guinea-Bissau are having super-refraction while all other part of the region
is dominated by sub-refraction. The area of sub-refraction was found to have
k-factor ranging between 1.24 and 1.32 (G value between -39.41 and -30.22
N-units/km) while the area of super-refractivity have k-factor of between 1.34
and 1.42 (G value of between -46.76 and -40.33 N-units/km). A very thin region
of normal refractivity (k = 1.33, G = -40 N-units/km) could be observed on
Figure 6B passing through southern Niger, part of north-western Nigeria,
Burkina Faso, Cote de Voire, Liberia, Sierra Leone, Guinea and south-western
Mali.

Figure 7: Dry Season map of (A) West African refractivity gradient G map and
(B) earth radius factor k map, for the first 1.5 km of the atmosphere.

By the dry season months of November, December, January, February, March and
April, the West African region is almost dominated by sub-refraction with only
small region of Senegal, Gambia, Guinea Bissau and Guinea showing some level
of super-refraction, there is also a very tiny region of normal refraction as
shown on Figure 7B, during the West African dry season the atmosphere usually
contains very little amount of water vapour. The region dominated with sub-
refraction was found to have k-factor within 1.22 and 1.32 (G value of between
-39.09 and -27.93 N-units/km) while the region of super-refraction was found
to have k-factor of between 1.34 and 1.42 (G value of between -46.20 and
-40.11 N-units/km).

Figure 8: Time series plot of refractivity gradient G and the k-factor for
four locations
($16~{}^{o}N,10~{}^{o}W;6~{}^{o}N,5~{}^{o}W;6~{}^{o}N,5~{}^{o}E;16~{}^{o}N,10~{}^{o}E$)
within the West African region.

A ten-year monthly average of k-factor and refractivity gradient G for four
locations in West Africa region is presented on Figure 8. The four locations
are chosen such that two are on either sides of the prime meridian
($5~{}^{o}E,10~{}^{o}E$ and $5~{}^{o}W,10~{}^{o}W$) as well as two at the
lower latitude and two at the higher latitude ($6~{}^{o}N,16~{}^{o}N$). It
could be observed on Figure 8 that the k-factor rises through the wet season
while the refractivity gradient falls through the wet season, following
variation in moisture through the seasons. Figure 8 confirms the fact that the
lower latitudes in West Africa generally have super-refraction while the
higher latitudes have sub-refraction.

##  4 Conclusion

Study of the radio refractivity, which is the change in the path of radio
signal as it passes through a medium, is an important part of modern
technology as it enhances modern communication, especially close to the
earth’s surface such as point-to-point terrestrial microwave communication,
terrestrial radio and television communication, wifi signals, and so on, also
important in the study of satellite signal propagation through the atmosphere.

The vertical variation of radio refractivity in the region was found to follow
an exponential model, that is, refractivity decreases exponentially with
height. The refractivity was also found to be higher at locations closer to
the ocean than locations closer to the desert. The current study has focused
on West African region within the latitudes $4~{}^{o}N$ and $18~{}^{o}N$, as
well as longitudes $18~{}^{o}W$ and $18~{}^{o}E$. The lower latitudes were
found to be dominated by super-refraction while the higher latitudes were
dominated by sub-refraction. The refractivity variation was found to be in
synchronization with the seasonal variation in atmospheric moisture across the
region.

## Acknowledgements

The authors wish to acknowledge the Radio Occultation Meteorology Satellite
Application Facility (ROM SAF) and Climate Monitoring Satellite Application
Facility (CM-SAF) of the EUMETSAT (data doi.org/10.5676/EUM_SAF
_CM/WVT_ATOVS/V001) for making their data freely available for use.

## References

  * [1] H. W. Ko, J. W. Sari & J. P. Skura, ”Anomalous microwave propagation through atmospheric ducts”, Johns Hopkins APL Technical Digest, 4(1) (1993) 12–26.
  * [2] ITU-R, ”The radio refractive index: its formula and refractivity data. ITU Radio Communication Assembly”, ITU-R P—Series, 453-11 (2015).
  * [3] S. M. Babin, G. S. Young & J. A. Carton, ”A new model of the oceanic evaporation duct”, Journal of Applied Meteorology, 36 (1997) 193–204.
  * [4] A. T. Adediji, M. O. Ajewole & S. E. Falodun, ”Distribution of radio refractivity gradient and effective earth radius factor (k-factor) over Akure, South Western Nigeria”, Journal of Atmospheric and Solar-Terrestrial Physics, 73 (2011) 12–26. doi: 10.1016/j.jastp.2011.06.017.
  * [5] M. P. M. Hall, ”Effect of the troposphere on radio communication, IEEE Electromagnetic Wave Series”, Peter Peregrinus Ltd. United Kingdom, 73 (1989) 105–116.
  * [6] T. J. Afullo, T. Motsoela & D. F. Molotsi, ”Refractivity gradient and k-factor in Botswana”, Radio Africa, 73 (1999) 107–110.
  * [7] A. T. Adediji & M. O. Ajewole, ”Vertical profile of radio refractivity gradient in Akure South-West Nigeria”, Progress in Electromagnetics Research C, 4 (2008) 157–168.
  * [8] O. A. Aremu, L. O. A. Oyinkanola, A. Akande & W. A. Azeez, ”Effects of radio refractivity gradient and k-factor on radio signal over Ibadan, south western, Nigeria”, Progress in Electromagnetics Research C, 5(6) (2018) 248–253.
  * [9] J. L. Davis, G. Elgered, A. E. Niell & C. E. Kuehn, ”Ground-based measurement of gradients in the wet radio refractivity of air”, Radio Science, 28(6) (1993) 1003–1018. doi: 0048-6604/93/93RS-01917508.00.
  * [10] O. A. Aremu, L. O. A. Oyinkanola, A. Akande & W. A. Azeez, ”A Tropical Model for Analyzing Radio Refractivity: Selected Locations in North Central, Nigeria”, Proceeding of the 2015 International Conference on Space Science and Communication (IconSpace), 10-12 August 2015, Langkawi, Malaysia (2015).
  * [11] O. L. Ojo, J. S. Ojo & P. Akinyemi, ”Characterization of secondary radioclimatic variables for microwave and millimeter wave link design in Nigeria”, Journal of Radio and Space Physics, 46 (2017) 83–90.
  * [12] O. O. Fashade, T. V. Omotosho, S. A. Akinwumi & K. P. Olorunyomi, ”Refractivity gradient of the first 1km of the troposphere for some selected stations in six geo-political zones in Nigeria”, IOP Conf. Ser.: Mater. Sci. Eng., 640 (2019) 83–90. doi: 10.1088/1757-899X/640/1/012087.