ing real structures on various scales and with differ-
ent shapes: there is a trade-off between localisation
accuracy and noise sensitivity. Nonlinear smoothing
has been developed to overcome these shortcomings,
which tend to preserve important features along with
noise removal during smoothing (Narendra, 1981;
Brownrigg, 1984; Perona and Malik, 1990; Nitzberg
and Shiota, 1992). In this context, the need for an
adaptive approach to cope with inhomogeneities in
images is well recognised (Saint-Marc and Medioni,
1991; Gomez, 2000; Grazzini et al., 2004). Strong
relations have been established between a number of
widely-used adaptive ﬁlters for digital image process-
ing (Mr
´
azek et al., 2006). The most common strategy
is to locally vary the kernel over image regions ac-
cording to their contents (Takeda et al., 2007): around
a pixel x to be updated, one uses a kernel K = K
x
with proper weights depending on the actual image
variability in the neighbourhood of x. A critical issue
is then how to measure image variability for gener-
ating the weights of the kernel. A possible approach
consists in adapting the local effect of the ﬁlter by
using both the location of the nearby samples and
their luminance values. That is, the proposed ker-
nel K
x
takes into account two factors: spatial dis-
tances |x −y| and tonal distances |f (x) − f (y)|. In-
troducing a tonal weight, the mixing of different lu-
minance ’populations’ is prevented. Such approach is
known as the bilateral ﬁltering technique introduced
by (Tomasi and Manduchi, 1998) as an intuitive gen-
eralisation of the Gaussian convolution. In fact, it
is clear from the literature that there is a direct re-
lationship between this technique and other nonlinear
smoothing methods (van den Boomgaard and van de
Weijer, 2003; Mr
´
azek et al., 2006): anisotropic dif-
fusion (Barash, 2002), local-mode ﬁnding (van den
Boomgaard and van de Weijer, 2002) or mean-shift
analysis (Comaniciu and Meer, 2002).
Following the line of thought of (Tomasi and Man-
duchi, 1998), this paper introduces a new algorithm
for edge-preserving smoothing of natural images as a
preprocessing stage in feature extraction and/or image
classiﬁcation. The present method exploits an image-
dependent approach derived from concepts known in
mathematical morphology (Soille, 2004). It is based
on the deﬁnition of appropriate geodesic masks and
the local estimation of pairwise geodesic time func-
tions within these masks (Lantu
´
ejoul and Maison-
neuve, 1984; Soille, 1994a; Soille, 1994b; Soille and
Grazzini, 2007). Likewise bilateral ﬁltering, the main
idea is to associate with each pixel a weighted con-
volution of sample points within an adaptive neigh-
bourhood, where the weights depend not only on the
spatial distances but also on the tonal distance to the
centre pixel. With respect to other nonlinear tech-
niques, which often involve iterative operations (Per-
ona and Malik, 1990; Comaniciu and Meer, 2002),
this approach presents the advantage of not depending
upon any termination time. Moreover, it enables to
determine adaptively, directly from the unsmoothed
input data, the neighbouring sample points and the
associated weights. Namely, by designing relevant
geodesic masks, the geodesic time functions com-
puted from a single pixel provide a similarity measure
of the twofold spatial and tonal information in the lo-
cal neighbourhood of this pixel.
The rest of the paper is organised as follows. In
the next section, we recall the fundamental notions
of geodesic path, geodesic mask and geodesic time
known in mathematical morphology. In section 3, we
introduce a new ﬁlter for image smoothing based on
the estimation of local geodesic time on the gradient
magnitude. In section 4, we propose an alternative al-
gorithm based on the calculation of the geodesic time
on image variations able to deal with the presence of
noise. In section 5, we present and discuss some re-
sults and also compare the approach with other ex-
isting techniques. A conclusion regarding the current
results and a description of future foreseen develop-
ments are presented in section 6.
2 GEODESIC TIME ON
GREYLEVEL IMAGES
Geodesic transforms are classical operators in image
analysis (Lantu
´
ejoul and Maisonneuve, 1984; Verwer
et al., 1989). The geodesic distance between two pix-
els of a connected set (typically, a binary image) is de-
ﬁned as the length of the shortest path(s) linking these
points and remaining in the set (the so-called geodesic
paths). This idea can been generalised to greylevel
images (images with single-valued luminance) using
the geodesic time on geodesic mask (Soille, 1994b;
Ikonen and Toivanen, 2007). In such case, the image
is then treated as a ’height map’, i.e. a surface embed-
ded in a 3D space, with the third coordinate given by
the greylevel values. The geodesic paths are then con-
strained to the surface of the height map: typically, the
path between two close pixels can be long, if there is
a high ’ridge’ in the greylevel map between them.
If we consider a greylevel image as an inte-
grable function g (designed as the geodesic mask), the
time τ
g
(P ) necessary for travelling on a path P de-
ﬁned on the domain of deﬁnition of g is expressed as
the integral of g along P (Soille, 1994b):
τ
g
(P ) =
Z
P
|g(s)|ds . (2)
EDGE-PRESERVING SMOOTHING OF NATURAL IMAGES BASED ON GEODESIC TIME FUNCTIONS
21