ImageAnalStereol2009;28:45-61
OriginalResearchPaper
LABELING OFN-DIMENSIONALIMAGESWITHCHOOSABLE
ADJACENCYOF THE PIXELS
KAISANDFORT1ANDJOACHIMOHSER2
1UniversityofKarlsruhe,IWRMM, Engesserstraße6,76131Ka rlsruhe,Germany;2UniversityofApplied
Sciences,Darmstadt,Sch¨ offerstraße3,46295Darmstadt, Germany
e-mail: sandfort@math.uni-karlsruhe.de,jo@h-da.de
(AcceptedFebruary18,2009)
ABSTRACT
Thelabelingofdiscretizedimagedataisoneofthemostesse ntialoperationsindigitalimageprocessing. The
notionsofanadjacencysystemofpixelsandthecomplementa rityoftwosuchsystemsarecrucialtoguarantee
consistency of any labeling routine. In to date’s publicati ons, this complementarity usually is defined using
discrete versions of the Jordan-Veblen curve theorem and th e Jordan-Brouwer surface theorem for 2D and
3D images, respectively. In contrast, we follow here an alte rnative concept, which relies on a consistency
relationfortheEulernumber. Thisrelationandallnecessa rydefinitionsareeasilystatedinauniformmanner
for then-dimensional case. For this, we present identification and c onvergence results for complementary
adjacencysystems,supplementedbyexamplesforthe3Dcase . Next,wedevelopapseudo-codefora general
labeling algorithm. The application of such an algorithm sh ould be assessed with regard to our preceding
considerations. Abenchmarkandathoroughdiscussionfinis hourarticle.
Keywords:adjacencyoflattice points,complementarity,c onnectivity,labeling,runlengthencoding.
INTRODUCTION
Labeling of connected components (‘objects’) is
one of the most important tools of image processing.
It is the basis for the generation of object features as
wellasofsomekindoffiltering, i.e.,removingofnoisy
objects or holes in objects. The criteria for removing
an object or a hole can be chosen extremely flexible
based on the object features. The task of labeling
(object filling, region detection) is to assign labels
(mostly unsigned integers) to the pixels in such a way
that all pixels belonging to a connected component
of the image are assigned the same label, and pixels
belongingto differentcomponentsgetdifferentlabels.
Due to its importance in image processing, there is
plenty of literature about labeling and techniques to
control (and improve) the processing and memory
demands. These can be tremendous for large images,
which frequently arise in practice. Once again, the
challenge lies in finding a compromise between
usability, flexibility,andefficiency.
The prototype of labeling algorithms is the simple
and well-known Rosenfeld-Pfaltz method (Rosenfeld
and Pfaltz, 1966; Klette and Rosenfeld, 2004). Here,
the image is scanned until a pixel xkis found that
has not yet been labeled. If there is no neighboring
pixel labeled w.r.t. the chosen adjacency system, a
new label is chosen for xk. Otherwise, if there is a
neighboring pixel xjwith the label ℓj, the label ℓj
is assigned also to xk. In the case of more than one
neighboring pixels which have different labels, thisis noted in a table of pairs of equivalent labels. A
connected component is finally identified as the set
of pixels belonging to the same equivalence class
of labels. The table of pairs can become very large
and may lead to problems in finding the equivalence
classes; the time complexity of a corresponding
algorithm is O(mlogm)withmbeing the number of
pairs. Thus, several variants of the Rosenfeld-Pfaltz
methodhavebeendevelopedwhichemploytechniques
to keepmassmallaspossible.
However, a reduction in complexity is achieved
also by other strategies. For example, in Dillencourt
et al.(1992) the authors start their analysis with
describing a combination of the Union-Find method
with weight-balancing and path-compression, which
yields O(˜mα(˜m)); here, αis the inverse of
Ackermann’s function and ˜ mdenotes the number of
FindandUnionoperationsandisrelatedtothevariable
mabove. On the one hand, the function αgrows
extremely slowly, and on the other hand, the quantity
˜mdoes not necessarily evolve on pixel level (as m)
since the Union-Find method can be equally applied
to other image representations with e.g., bintrees or
run lengths (instead of pixels) as modules. The same
holdsforextensionsof Union-Find (Dillencourt etal.,
1992)andrelatedalgorithms.Thefactthatsomeimage
representations,inparticularrunlengthencodings,are
compatible with and facilitate the labeling w.r.t. an
adjacencysystemprovidesthebasisforouralgorithm.
Given the input image as a pixel array, a conversion
into a different representation might be done as a
45
SANDFORT KET AL:Labelingofn-dimensionalimages
preprocessing in a labeling scheme. Its purpose could
be a decomposition (Aguilera et al., 2002) or a more
compact representation (our method) of the image to
allowanefficientdataaccessand/ortoreducememory
demands.
Many classical labeling methods, including the
one by Rosenfeldand Pfaltz, are two-pass-techniques.
Here, in the first pass preliminary labels are assigned
and equivalencesbetween labels are registered, and in
the second pass these correspondences are resolved
into equivalence classes. Either representatives of
these classes ( e.g., their smallest elements) or
consecutive integers assigned to them are then set
as the final labels. The resolving step is a critical
issue, and several methods have been proposed to
handleit.AnefficientsolutionisexplainedinThurfjell
et al.(1992). However, the basic idea to maintain
and merge equivalences in a 1D array originated in
a work by Hoshen and Koppelman (1976). Hoshen
(1998) later on examined its use in image processing.
An alternative and general solution, with special
consideration of dismissing and reusing preliminary
labels, is given in Dillencourt et al.(1992) (integrated
in aUnion-Find method), along with a detailed proof
of correctness. While in the latter the final label is
determined only in a second run, the methods cited
before can be used to merge label classes on-the-fly,
i.e., as soon as an equivalence is detected, and by that
allow a direct mapping of a preliminary label to the
correspondingfinalonerightafterthefirstrunthrough
the inputimage.
A very different concept is followed by recursive
labeling methods. Such methods completely avoid
the setup of a mapping between preliminary
and final labels and identify a complete object
(connectedcomponent)atonce.Theinformationabout
connectedness is usually carried in the algorithms by
design, rather than provided as a parameter or data
structure. An early account of recursive detection
of connected components is given in Hopcroft and
Tarjan (1973). The underlying principle resembles
that of filling algorithms. Recursive techniques can
achieve an excellent performance as they often are
accessible to the optimization potential of a compiler
andaCPU.However,sincetherecursiondepthheavily
depends on the contents of the input image and so
can not be predicted, they have to tackle the risk of a
stack overflow. A smart approach to this is described
by Mart´ ın-Herrero in Mart´ ın-Herrero (2004), see
also Mart´ ın-Herrero (2007). It relies on a quasi-run
lengthencodingin onedimensionandrecursioninthe
remainingdimensionsofthe inputimage,whichgains
asignificantdecreaseoftherecursiondepthcompared
torecursioninalldimensions.Theresultingalgorithmyields performance values which seem to be (among)
the best of today’s labeling algorithms and will serve
asabenchmarkin ourpaper.
We now explain what we mean by quasi-run
length encoding above and put the bridge to our
labeling method.We assumea binary inputimage and
consider the labeling of its foreground pixels. In the
code by Mart´ ın-Herrero, starting with a single pixel
P, all foreground pixels which are consecutive in a
given direction and trivially connected to Pare found
by iterative scanning and assigned a label differing
from the foreground value. Afterwards, the same is
recursively performed for all foreground neighbors
(accordingtothenotionofconnectedness)ofallpixels
oftherun.Eachneighboractsasthestartingpixel Pof
thenextstageoftherecursion.Inthisway,allpixelsof
a connectedcomponentare reached and every pixel is
assigned a label only once. However, in searching for
the foreground neighbors of all members of a run, a
lot of pixels are accessed more than once. Although
the query of a label is a simple and fast operation,
this could be seen as a (small) shortcoming. It is not
possibleheretorestrictthecheckoftheneighborstoa
fewofthepixelsofarunsinceotherwisenotallpixels
of an object might be captured by recursion. Anyway,
this hybrid method is a single-pass-technique in the
sensethat the input image as a whole has to be passed
only once in order to labelc orrectly all connected
components.
For our algorithm, we take a perspective which is
somehow complementary to the one just described.
We start with an iterative run length encoding
of all foreground pixels in the input image. This
preprocessing has been mentioned before as a
conversion from one image representation into
another. It guarantees that all foreground pixels are
registered a priori and makes it possible to detect
all correspondences by checking only the boundaries,
i.e., the start and end point of each run. The overall
labeling algorithm permits a proper control of the
stack and memory demands. We also believe that
the representation of the image as an array of run
lengthsisadvantageousforafurtherimageprocessing
as well as analysis, e.g.a fast extraction of object
features. On the flipside, however, it lacks some of
the parallel processing and real time features of the
method by Mart´ ın-Herrero if the image is given as a
pixel array. Both methods exploit a property of any
reasonable definition of adjacency, namely that pixels
which are consecutive in a direction are connected.
This is reflected below in the inclusion Eq. 3. Here,
all admissible directions are implicitly defined by the
lattice orthepixelarray, respectively.
46
ImageAnalStereol2009;28:45-61
Finally, we remark that the introduction in
Thurfjell et al.(1992) contains a short survey of
different labeling algorithms, and Chapter 6 in Ritter
and Wilson (2001) formalizes labeling operations on
the basisofimagealgebra.
The outline of the paper is as follows. First we
give a short introduction to lattices (point lattices,
grids). For a general definition of pixel neighborhood,
we follow the approaches of Ohser et al.(2002;
2003); Schladitz et al.(2006) and state the concept
of adjacency on n-dimensional lattices and pairs of
complementary adjacency systems in. An important
insight is that, in 3D, complementarity according to
a consistency relation for the Euler number differs
from the complementarity as conventionally defined
bythediscreteJordan-Brouwertheorem.Furthermore,
werecallthenotionofconnectednessinthecontinuous
case and give a definition of connected components
in a lattice (discrete case) w.r.t. to a given adjacency
system. This finishes our theoretical investigation.
Beforeexplainingthelabelingalgorithm,weintroduce
somebasicmathematicalobjectstohandleimagesand
intendtoclarifythelinkbetweenthebackgroundfrom
the previous sections and a practical implementation.
Afterwards, we explain a run length encoding and
a labeling procedure by means of pseudo-codes. For
a convenient illustration, important and frequently
used functionality is encapsulated here in auxiliary
routines. A further section is devoted to a benchmark
and a test of our algorithm on application data from
material testing. We finally address some applications
and possible extensions of the algorithm in a short
conclusion.
ADJACENCYOFHOMOGENEOUS
LATTICESANDEULERNUMBER
Ann-dimensional homogeneouslattice is a subset
Lnofthen-dimensionalEuclideanspace Rnwith
Ln={x∈Rn:x=n
∑
i=1λiui,λi∈Z}=UZn,(1)
whereu1,...,un∈Rnform a basis of Rn,U=
(u1,...,un)is the matrix of column vectors, and Zist
the set of integers. By C=U·[0,1]nwe denote the
unit cell of Ln. LetF0be the set of vertices of a
polyhedron,in particular F0(C) =U·{0,1}n.
In the literature, the adjacency of lattice points
(or pixels) is usually characterizedby a neighborhood
graphΓ, and the complementarity of adjacencies
is defined via the Jordan-Brouwer surface theorem
(Lachaud and Montanvert, 2000). Here we use analternative concept of adjacency systems based on
the Euler number of a discretization, thoroughly
introduced in Nagel et al.(2000) and Ohser et al.
(2002;2003).
DISCRETIZATIONWITHRESPECT TO
AN ADJACENCY SYSTEM
LetLnbe a lattice with the basis {u1,...,un}and
the unit cell C. The vertices of Care indexed, and
we write xj=∑n
i=1λiui,λi∈ {0,1}, with the index
j=∑n
i=12i−1λi.Clearly,theunitcell Chas2nvertices,
xi∈F0(C),i=0,...,2n−1. In a similar way we
introduce the index of a subset ξ⊆F0(C). Let 1
denotetheindicatorfunctionofaset, i.e.1(x∈ξ) =1
ifx∈ξand 1(x∈ξ) =0 otherwise. The index ℓis
assigned,andwewrite ξℓif
ℓ=2n−1
∑
j=02j·1(xj∈ξ), (2)
for allξ⊆F0(C), where ℓtakes the integer values
between 0 and ν=22n−1. Theξℓare local pixel
configurationoftheforegroundofabinaryimage,and
ν+1is thenumberofdifferentconfigurations.
We introduce the convex hulls Fℓ=convξℓ
forming convex polytopes with Fℓ⊆CandF0(Fℓ)⊆
F0(C),ℓ=1,...,ν. LetFj(F)denote the set of all
j-dimensional faces of a convex polytope F. For a set
Fof convex polytopes we set Fj(F) =/uniontext{Fj(F):
F∈F}. Now we are able to equip the lattice Ln
with a (homogeneous) adjacency system defining the
neighborhoodoflattice points.
Definition 1 LetF0⊆ {F0,...,Fν}be a set of convex
polytopesF ℓ=convξℓ, and let Fbethe union overall
lattice translationsof F0,thatis F=/uniontext
x∈LnF0+x. If
(i)/ 0∈F0,C∈F0,
(ii) ifF ∈F0thenFi(F)⊂F0fori=0,...,dimF,
(iii)if F i,Fj∈F0andconv(Fi∪Fj)/∈F0then Fi∩Fj,
Fi\Fj,Fj\Fi∈F0(where the bar denotes the
topologicalclosure),
(iv)ifF i1,...,Fim∈F0andF =/uniontextm
j=1Fijisconvexthen
F∈F0, m=2,...,ν,
then the system F0is called a local adjacency system
andFis said to be an adjacency system of the
lattice Ln.
ThepairΓ= (F0(F),F1(F))istheneighborhood
graphofFconsisting of the set F0(F)of nodes and
the set F1(F)of edges. The order of the nodes is
calledthe connectivity ofLn.
47
SANDFORT KET AL:Labelingofn-dimensionalimages
Inthesimplestcase,wheretheadjacencysystemis
generated from the unit cell C, the order of the nodes
is 2n, and we write F2n=/uniontext
x∈Ln/uniontextn
j=0Fj(C+x). The
maximum adjacency system consisting of the convex
hullsofallpointconfigurationsprovidesa κ-adjacency
withκ=3n−1,Fκ=/uniontext
x∈Ln{F0+x,...,Fν+x}.
Notice that for all adjacency systems FonLnthe
inclusion
F2n⊆F⊆Fκ (3)
holds. Now we recall the adjacency systems on L3
consideredindetailin Ohser etal.(2002;2003).
–6-adjacency. The6-adjacencyisusedasastandard
in image processing. It is generated from the unit
cellC,F6=/uniontext
x∈L3/uniontext3
j=0Fj(C+x).
–14.1-adjacency. This adjacency system is
generated from the tessellation of Cinto the 6
tetrahedra F139,F141,F163,F177,F197, andF209
beingtheconvexhulls oftheconfigurations
,,,,,.
Thatis F0, consistsofall j-faces ofthetetrahedra,
j=0,...,3,andtheirconvexunions.Theedgesof
the corresponding neighborhood graph Γare the
edges ofC, the face diagonals of Ccontaining the
origin 0,the spacediagonalof Ccontaining0,and
alltheirlattice translations.Theorderofthenodes
ofΓis14.
–14.2-adjacency. The 14.2-adjacency system is
generatedfrom the tetrahedra F43,F141,F147,F169,
F177, andF212,whichare theconvexhulls of
,,,,,.
The corresponding neighborhood graph Γdiffers
from that one for 14.1 in the choice of one face
diagonalof Csuchthatit doesnotcontain0.
–26-adjacency. The system F26is the maximum
adjacenysystemon L3,F26=Fκ.
Based on the definition of adjacencywe introduce
theF-discretizationofaset.
Definition 2 The discretization X ⊓Fof a compact
subset X ⊂Rnw.r.t. a given adjacency system Fis
definedas the unionof all j-faces of the elementsof F
for whichallthe verticeshitX,i.e.,
X⊓F=/uniondisplay
{F∈F:F0(F)⊆X}.(4)
It is important to realize that in particular for
higher-dimensional images the connectivity of the
pixels and, hence, the labeling can heavily depend onthe choice of the adjacency system. The number of
neighbors of a pixel in an n-dimensional image can
rangefrom2 nto3n−1.Forthe2dcasetheconectivity
ranges from 4 to 8 while in the 3d case the extremal
choices are 6- and 26-connectivity. As a consequence
of the wide range of the number of neighbors,
the neighborhood of pixels should be chosen very
carefully in dependence of the dimensionality of the
image, the lateral resolution, the image data, and the
aims of processing and analysis (Schladitz et al.,
2006).
EULER NUMBER AND
COMPLEMENTARITY OFADJACENCY
SYSTEMS
Since the set X⊓Fforms a (not necessarily
convex) polyhedron, the number ♯Fj(X⊓F)of
elements of Fj(X⊓F)is finite and, therefore, the
Eulernumber χ(X⊓F)canbecomputedviatheEuler-
Poincar´ eformula,
χ(X⊓F) =n
∑
j=0(−1)j♯Fj(X⊓F).(5)
A local version of Eq. 5 is given in Schladitz et al.
(2006).
Itiswell-knownfrom imageprocessingthatifone
choosesanadjacencysystem Fonthediscretizationof
X, then there is implicitly chosen a system Fcon the
discretization of the complementary set Xc. In other
words, if the ‘foreground’ X∩Lnis connected w.r.t.
F, then the ‘background’ Xc∩Lnmust be connected
w.r.t. Fc. Forn>2 it is not sufficient to consider
connectivity, and further criteria have to be regarded.
In the following we introduce ‘complementarity’ by
meansoftheEulernumberofthediscretization X⊓F.
Noticethatthecomplement XcofXisunbounded.
Thus, the Euler number χ/parenleftbig
Xc/parenrightbig
of its topological
closure may be defined by Hadwiger’s recursive
definition(Schneider,1993,p.175).
Definition3 The pair (F,Fc)is called a pair of
complementary adjacency systems if (X⊓F)∩(Xc⊓
Fc) =/ 0and
χ(X⊓F) = (−1)n+1χ(Xc⊓Fc)(6)
holds for all compact X ⊂Rn. An adjacency
system Fis called self-complementary if χ(X⊓F) =
(−1)n+1χ(Xc⊓F)holdsfor allcompactX.
48
ImageAnalStereol2009;28:45-61
Forn=3 there exist 3 pairs of complementary
adjacency systems: (F6,F26),(F14.1,F14.1)and
(F14.2,F14.2). That is, the 6-adjacency is
complementary to the 26-adjacency and there are
known two self-complementary adjacency systems,
the 14.1-adjacency and the 14.2-adjacency, see Ohser
etal.(2002;2003).
Consider a lattice L3equipped with a
neighborhood graph Γ′= (L3,F1)where the system
of edges F1may consist of all edges and face
diagonals of the cells of L3. The order of the nodes
ofΓ′is 18 and, hence, the adjacency is called
the 18-adjacency (which is widely used in image
processing). The 18-adjacency is ‘Jordan-Brouwer-
complementary’ to the 6-adjacency generated solely
from the edges of the lattice cells (Lachaud and
Montanvert, 2000). However, Eq. 6 does not hold for
the pair (F18,F6)and hence it is not complementary
in the sense of Definition 3 (Schladitz et al., 2006).
This shows that in higher dimensions ( n>2) the
‘Jordan-Brouwer-complementarity’ differs from that
of Definition 3. The criteria in Definition 3 seam
to be stronger than that of the ‘Jordan-Brouwer-
complementarity’.
MULTIGRIDCONVERGENCE
Now we consider the relationship between the
Euler number of a compact set X⊂Rnand the Euler
number of its discretization. It can not be expected
thatχ(X) =χ(X⊓F)for all compact sets X, but if
Xhas a sufficiently smooth surface, the Euler number
ofX⊓Fconverges to the Euler number of Xfor
increasing lateral resolution. Here ‘smooth’ is defined
by morphologicalopeningandmorphologicalclosure.
The setXis called morphologically open w.r.t. a set
A⊂RnifXisinvariantw.r.t.openingwith A,X◦A=
X.Heretheopeningisdefinedby X◦A= (X⊖ˇA)⊕A,
andX⊖A= (Xc⊕A)cis the Minkowski subtraction.
Analogously, the set Xis called morphologically
closedw.r.t. AifX•A=XwhereX•A= (X⊕ˇA)⊖A
is the morphological closure with A. Morphological
regularity of Xmeans that there is an ε>0 such that
Xis morphologicallyopenaswell as morphologically
closedw.r.t. a ball Bεofradius ε.
Theorem1 Let(F,Fc)be a pair of complementary
adjacency systems on Ln. If X ⊂Rnis
morphologically closed w.r.t. all edges F ∈F1(F)
andmorphologicallyopenw.r.t.allF ∈F1(Fc),then
χ(X) =χ(X⊓F)andχ(Xc) =χ(Xc⊓Fc).
AproofisgiveninOhser etal.(2002).Noticethataset
Xfulfillingthelastconditionispolyconvexand,hence,its Euler number exists.However, this condition for X
is very strong, it depends on Fand, hence, it will not
be fulfilled in most applications. Thus, we consider a
more natural condition for X. LetBεbe a (small) ball
of radius ε. From Theorem 1 we obtain the following
lemma.
Lemma1 Let(F,Fc)be a pair of complementary
adjacency systems on Ln. Then a Fis an adjacency
systemon a Ln, a>0,andit is
lim
a→0χ(X⊓aF) =χ(X) (7)
for allcompactandmorphologicallyregularsetsX.
ThismeansthattheEulernumberisconvergentfor
morphologically regular sets (multigrid convergence).
Theproofofthislemmafollowsfromthefactthatif X
is morphologically regular, there exists an a>0 such
thatχ(X•F) =χ(X)for allF∈aFandχ(X◦F) =
χ(X)forF∈aFc, and choose an a>0 such that
F⊂Bεfor allF∈aF∪aFc.
CONNECTEDNESS
In order to describe a labeling algorithm, it is
necessary to introduce the notions of ‘connectivity’
and ‘connected component’. These are provided
by topology. Azriel Rosenfeld introduced a digital
topology on L2(Rosenfeld, 1970). He defined
connectednesson lattices andstateda discreteJordan-
Veblencurvetheorem.Thedefinitionofconnectedness
can simply be extended to n-dimensional lattices
(Lachaud and Montanvert, 2000). Here we introduce
connectednessbasedonadjacencyoflattice points.
CONTINUOUSCASE
Firstly we consider the continuous case and
introduceconnectivityfortheEuclideanspace Rn.The
connectedcomponentsofaboundedset X⊂Rncanbe
consideredastheequivalenceclassesof X⊆Rnw.r.t.
anappropriatelychosenequivalencerelation ∼defined
forpointpairsin Rn.
Definition 4 A set X ⊂Rnis said to be connected if
for all subsets X 1,X2⊆X with X 1∪X2=X it follows
thatX1∩X2/\e}atio\slash=/ 0orX1∩X2/\e}atio\slash=/ 0.
This definition of connectivity is closely related
to path-connectivity. A path in Rnis a continuous
mapping f:[0,1]/mapsto→Rn. Iff(0) =xandf(1) =y,
x,y∈Rn, thenfis calledapathfrom xtoy.
49
SANDFORT KET AL:Labelingofn-dimensionalimages
Definition 5 A non-empty set X is called path-
connected if for every x ,y∈X there exists a path f
fromx to ysuchthat f (·)⊆X.
It is well-known that every path-connected set
Xis also connected. Furthermore, if Xis open and
connected, it is also path-connected. Obviously, a
connected set Xis not necessarily path-connected.
For example, the curve of the function sin (1/t)is
connected,butnotpath-connected.Moreprecisely,the
setX={(t,sin1
t):t∈R\{0}}∪{(0,t):t∈[−1,1]}
is connected,butnotpath-connected(Schmitt, 1998).
We write x∼yfor path-connected points x,y∈
Rn. It can be shown that the binary relation ∼is an
equivalence relation, i.e.∼is reflexive, symmetric,
and transitive. The equivalence classes X1,...,Xmof
Xunder ∼are called path componentsof X.Formore
detailssee, e.g.,Armstrong(1997)andRotman(1993).
DISCRETECASE
Connectednessinadiscretizationiscloselyrelated
to adjacency of lattice points. Hence, we consider
a homogeneous lattice Lnequipped with a pair of
complementary adjacency systems (F,Fc). Letxand
ybe lattice points, x,y∈Ln. A discrete path from xto
yw.r.t.theadjacencysystem Fisasequenceoflattice
points (xi)m
i=0⊂Ln,m∈N, withx0=x,xm=y, and
[xi−1,xi]∈F,i=1,...,m.
A non-empty discrete set Y⊆Lnis called path-
connected w.r.t. Fif♯Y=1 or if for all pairs (x,y)∈
Y2withx/\e}atio\slash=ythereexistsadiscretepathw.r.t. Ffrom
xtoy. Connectedness w.r.t. FinYis an equivalence
relation.
Definition 6 LetLnbe a homogeneous lattice
equipped with an adjacency system F, and letY ⊆Ln
be a discreteset. The equivalenceclassesY 1,...,Ym⊆
Y, m≥1, defined through the connectedness w.r.t. F
arecalledthe connectedcomponentsofY.
We will use the notation YF={Y1,...,Ym}for the set
of equivalence classes of Yw.r.t. F. The following
Lemma links the equivalence classes of a set X⊂Rn
to theequivalenceclassesofits digitisation X∩Ln.
Lemma2 Let(F,Fc)be a pair of complementary
adjacencysystems on Ln, and let X be a compactand
morphologically regular subset of Rnwith the set of
equivalenceclasses {X1,...,Xm}under ∼. Then there
is aconstantb >0suchthat
(X∩aLn)aF={X1∩aLn,...,Xm∩aLn}(8)
for allawith 0<a<b.Proof.IfXis morphologicallyclosedw.r.t. Bε,ε>0,
thenXi∩Xj=/ 0 implies that
inf{/bardblxi−xj/bardbl:xi∈Xi,xj∈Xj}>ε.
Nowwechoose bsuchthatb
2(C⊕ˇC)⊂Bε,whereCis
theunitcellof Ln.ThenXi⊓aFandXj⊓aFaredisjoint
and,thus, there is no discrete pathw.r.t. Fconnecting
Xi∩aLnandXj∩aLn.
On the other hand, since Xiis morphologically
openw.r.t. Bε, then for eachpath finXiexistsa path
gwithf⊂g⊕Bε⊆Xi. Since (g⊕Bε)∩aLnis path-
connected w.r.t. aF, it follows that Xi∩aLnis path-
connectedw.r.t. aF, too. /square
From Lemma 2 it follows that for sufficiently
high lateral lattice resolution the equivalence classes
of(X∩aLn)are independent of the choice of the
adjacency system. However, this holds only for sets
Xwith sufficiently smooth surface. In general, the
equivalenceclassesof X∩Lndependon (F,Fc).
LABELINGWITH CHOOSABLE
ADJACENCY
Inthissection,weconsideranimplementationofa
generalandcustomizablelabelingalgorithm.Essential
for this and its performance is the exploitation of the
fact that runs are naturally connected and as such
part of the same component, cf.the inclusion Eq. 3.
In particular, this complies with all local adjacency
systemsgivenaboveforthecase n=3.Webeginwith
a short description of the labeling procedure. Then
we explain the necessaryvariables and data structures
andpresenttheproceduresaseasy-to-translateC-style
pseudo-codes.
BASICS
We anticipate that an image I⊂Lnis a
finite discrete set of lattice points (pixels) with
an associated binary-valued color mapping v:I/mapsto→
{BGCOL,FGCOL}, where BG COL≤0 stands for
the value of a background color, and FG COL is the
value of a foreground color. We simply speak about
foreground and background pixels (of I) and by that
meanthose elements xofIfor which v(x) =FGCOL
andv(x) =BGCOL, respectively. Let Y={x∈I:
v(x) =FGCOL}denotethesetofforegroundpixels.
Preceding the labeling, our algorithm does a
run length encoding of the foreground pixels of I
to identify Ybeforehand and to allow a compact
representation. As a runwe consider a set of
50
ImageAnalStereol2009;28:45-61
consecutive pixels in a certain direction of I. All
admissible directionsare indicated hereby the vectors
u1,...,unof the underlying lattice Ln, see Eq. 1. The
runlengthencodingassignseitherasingleplaceholder
label to all runs or a unique preliminary label to
each run. This depends on whether, in a subsequent
labeling, labels should be propagated to adjacent
runs or not. During the labeling, these labels are
systematically overwritten in a way which allows to
findthe connectedcomponentsin Y.
Using the notation introduced above, in formal
terms a labeling w.r.t. an adjacency system Fis a
mapping LF:I/mapsto→N0definedby
LF(x) =/braceleftBigg
f(i)forx∈Yi,
0 otherwise, i.e.x∈I\Y,
wheref:{1,...,m} /mapsto→Nis a fixed function with
f(i)/\e}atio\slash=f(j)fori/\e}atio\slash=j. For simplicity, we let f=idin
the following.
To ‘cluster’ the runs which belong to the same
connected component, equivalences between their
preliminary labels with regard to F(also called
‘correspondences’)havetoberegisteredandprocessed
such that finally all runs making up the i-th connected
component are assigned the unique integer f(i) =
ias a final label. The central factor in finding all
correspondences in an efficient way is the relation of
Eq. 3, which implies that the foreground neighbors
of only the start and endpoint of each run have to be
examined.
Now, we consider some basic mathematical
objects to handleimages.Let X⊂Rnbe as in Lemma
2,M⊆¯M={1,...,n}andUMbe a matrix consisting
of the column vectors ujwithj∈M. In particular, it
holdsUZn=U¯MZn. For a vector v, letvjdenote its
j-th element.Wedefine
the window:W =URnwithRn={r∈Rn: 0≤r(j)≤
d(j)for1≤j≤n}forafixed d∈Nn,
theimage: I=Ln∩W=U(Zn∩Rn)where Ln=UZn,
thecolormapping:v :I/mapsto→ {BGCOL,FGCOL}with
v(x) =/braceleftBigg
FGCOL for x∈I∩X
BGCOL for x∈I\X,
the setofforegroundpixels:Y =I∩X,
projections:P M(Ln) =UMZ♯M⊆Ln, and
a subimage:S M,x= (x+PM(Ln))∩W⊆Iforx∈I.Here, ♯Mdenotes the number of elements in M.
Associatedwithanyset A⊆Iisthecolormapping v|A.
Asspecialcases,wementionthata sliceisasubimage
SM,xwith♯M=2, and alineis a subimage SM,xwith
♯M=1. Clearly, every line is the intersection of two
slicesSM1,xandSM2,xwith♯(M1∩M2) =1. Note that
SM,x=SM,yforx−y∈PM(Ln).
Thefollowing aspectswill beimportant:
–For an-dimensional image Ithere are 2nn!
possibilities to scan through Ialong its nlattice
directions. The factorial in this term originates
from the scanning order of the direction indices,
and the factor 2naccounts for the orientations of
passing the directions. From now on, we assume
that these orientations coincide with those of the
basis vectors u1,...,unsuch that n! possibilities
remain. The scanning order is given by the ranks
ofthedirections,wherethefirstscanningdirection
hasrank 0,thesecondonerank 1, etc.
–Runs of pixels are detected and coded in the rank
0-direction. For an anisotropic set X∩W(with
discrete analogon Y), the processing speed of our
labeling algorithm depends on the choice of the
rank 0-direction since it determines the number of
runs.
Now we describe the main data structures which
weuseinourpseudo-codesfortherunlengthencoding
and the labeling. To simplify understanding the code,
scalar variables begin with a small letter and data
structures including vectors beginwith a capital letter.
All vectors are assumed to provide the methods
add()toaddanelement, at()toaccesstheelement
with the index given by the argument (starting from
0),init() to set all elements to the value of the
argument, resize() toresizethevectortothelength
given by the argument and size() to query the
number of elements. The i+1-th element can be
accessed alternatively by the [i]-operator, following
the nameofthe vector.
The main vectors, arrays and structures used in
our run length encoding and labeling procedures
areImage,ImgSize ,Rank,Index,RleRun,
RleLine andNeighbors .
ByImagewedenotethepixelarrayfortheimage
data{(x,v(x)):x∈I}, precisely the entry for the
pixelx∈Ihas the value v(x). We assume that the
entry is accessed by Image[x1]...[xn]where here
(x1,...,xn)is the coordinate vector for xw.r.t. the
basis {u1,...,un},i.e. x =x1u1+...+xnun. In an
implementation of an algorithm capable to handle
images of an arbitrary dimensionality nwhich is
unknown initially, the addressing of pixels is a bit
51
SANDFORT KET AL:Labelingofn-dimensionalimages
more complex since the above style is inappropriate
for unknown n. However, this is not an issue of
interest here, so for illustration we choose the simple
form. The vector ImgSize stores the dimensions
ofImage. It is the pendant of the vector din the
definition of the window from the model. Rankis a
vector of length nand stores the ranks of the lattice
directions.Thevalue Rank[i] withi∈ {0,...,n−1}
is the rank of the direction with index i(indicated by
ui).Thescanningorderis Rank[0] ,...,Rank[n-1] .
Indexis the complementary vector (of length n)
ofRank, which maps the rank of a direction
to its index, i.e.Index[Rank[i]] = i and
Rank[Index[r]] = r for 0 ≤i,r<n. The
structure RleRun, holding the data of a single run,
consists of the attribute labeland the two-element
vectorPos. The value Pos[0] is the coordinate of
the start pixel of the run in the rank 0-direction, and
Pos[1] is the corresponding coordinate of its end
pixel.Theprocedurefortherunlengthencodingwrites
a preliminary label into label. Finally, RleLine
is a vector of RleRuns and represents a run length
encodedline S{i0},xoftheimage I.Here,i0denotesthe
index of the rank 0-direction. The size of RleLine
depends on the contents of the image, the index i0
and the line with representative x∈I. The array
Neighbors representsavectorofrelativecoordinate
vectors for the neighbors of a pixel according to
the governing notion of connectedness. This means
that the neighbors of a pixel x∈Ihave coordinate
vectors of the form (x1,...,xn)+Neighbors[j] ,
wherej∈ {0,...,κ−1},κisthenumberofneighbors
andNeighbors[j] is a coordinate vector of length
nwithNeighbors[j][k] ∈ {−1,0,1}fork∈
{0,...,n−1}. The order in which the neighbors are
stored in Neighbors does not matter. This array
represents the neighborhood graph of an adjacency
system accordinge to Definition 1. Due to the Eq. 3,
Neighbors is expected to contain the 2 nvectors
(0,...,vk,...,0)withk∈ {1,...,n}andvk∈ {−1,1}beingthe k-thelement.
THERUN LENGTH ENCODING
ALGORITHM
Besides the data structures from above, the
implementation of our run length encoding needs the
following. Entities marked as ‘global’ are available in
allroutines.Thevariables n,FG_COL andBG_COL as
wellasthevectors RankandIndexareglobal.
IND– Boolean variable which indicates whether a
run of foreground pixels is active ( true) or not
(false)
LABEL_PROP [global] – Boolean variable which
indicateswhetherlabelpropagationis switchedon
(true)oroff( false)
totNumLin –the totalnumberoflines
prelLabel – variable for assigning a
unique preliminary label to each run if
LABEL_PROP == false
CurPixCs – vector holding the coordinates
of the current pixel (w.r.t. the original
order of directions), i.e.the current pixel is
Image[CurPixCs[n-1]]..[CurPixCs[0]] .
LinePos – vector storing the position of the current
line in the image w.r.t. the scanning order, i.e.
CurPixCs[Index[r]] = LinePos[r-1]
for0<r<n
RleData [global] – a n−1-dimensional
array of RleLine s with dimensions
(ImgSize[Index[n-1]] ,...,
ImgSize[Index[1]] ), stores the run length
encodedversionof Image
We are now prepared to turn to the pseudo-code
for the procedure DoEncoding() , which performs
therun lengthencoding.
void DoEncoding()
{
bool IND;
int i, r, totNumLin, prelLabel;
vector<int> CurPixCs(n), LinePos(n-1);
RleRun Run;
RleLine *pLine;
// set variables
CurPixCs.init(0);
LinePos.init(0);
totNumLin = 1;
Index[Rank[0]] = 0;
FOR i FROM 1 TO n-1 DO
{
Index[Rank[i]] = i;
52
ImageAnalStereol2009;28:45-61
totNumLin *= ImgSize[i];
}
// allocate memory for the array RleData
// The size of the i-th dimension (i = 1,...,n-1) of RleData is the
// image size in the direction with rank n-i.
RleData = new RleLine[ImgSize[Index[n-1]]]..[ImgSize[I ndex[1]]];
// main loop
DO
{
pLine = &RleData[LinePos[n-2]]..[LinePos[0]];
// set the coordinates of the current pixel
FOR r FROM 1 TO n-1 DO
CurPixCs[Index[r]] = LinePos[r-1];
CurPixCs[Index[0]] = 0;
// check whether a run starts at the first pixel in the current line
IF Image[CurPixCs[n-1]]..[CurPixCs[0]] == FG_COL THEN
{
IND = true;
Run.Pos[0] = 0;
}
ELSE
IND = false;
FOR i FROM 1 TO ImgSize[Index[0]]-1 DO
{
CurPixCs[Index[0]]++;
IF Image[CurPixCs[n-1]]..[CurPixCs[0]] == FG_COL THEN
{
// if IND == true, then current run continues, nothing to do
// a new run starts
IF !IND THEN
{
IND = true;
Run.Pos[0] = i;
}
}
ELSE
{
// current run ends
IF IND THEN
{
IND = false;
Run.Pos[1] = i-1;
prelLabel++;
Run.label = LABEL_PROP? -999:prelLabel;
pLine->add(Run);
}
// if IND == false, then BG_COL-run continues, nothing to do
}
}
// the end of the current line in the rank 0-direction is reach ed
// check whether some run is active and end it
IF Image[CurPixCs[n-1]]..[CurPixCs[0]] == FG_COL THEN
{
53
SANDFORT KET AL:Labelingofn-dimensionalimages
IF !IND THEN
Run.Pos[0] = ImgSize[Index[0]]-1;
Run.Pos[1] = ImgSize[Index[0]]-1;
prelLabel++;
Run.label = LABEL_PROP? -999:prelLabel;
pLine->add(Run);
}
// update the line position, according to the scanning order
r = 0;
WHILE r < n-1 AND ++LinePos[r] >= ImgSize[Index[r+1]] DO
LinePos[r++] = 0;
}
WHILE LinePos[n-2] < ImgSize[Index[n-1]]
}
THELABELING ALGORITHM
Next, we consider the actual labeling of the
image data, which is based on a preceding run
length encoding. We have already stated that all
correspondences between preliminary labels can be
found by testing only the start and the endpoint of
runs with their respective neighbors according to the
chosenadjacency.Thenumberofcorrespondencescan
be kept very small if not a unique preliminary label is
assignedtoeachrunbeforehand(duringtherunlength
encoding), but instead labels are propagated to all
adjacent,non-labeledruns(havingstilltheplaceholder
label -999), and a new preliminary label is given only
to a run which has not already received one, while
scanning the image in the specified order. If labels
are propagated, only correspondences between a new
preliminary label and a previously propagated one
haveto be registeredandresolvedlater on.Otherwise,
the number of preliminary labels and the number of
arising correspondences usually is much higher. A
recursive labeling approach like that from Mart´ ın-
Herrero (2004) is based on an extremal form of label
propagation in the sense that the propagation may
emerge from every point of a run and is started
again in every line which received a propagated label.
In our algorithm, it is restricted to the start and
endpointofarunandstoppedattheadjacentruns.The
number of preliminary labels determines the size of
the vector LabelMap , whose entry indices represent
the preliminary labels. This vector is processed by
means of the correspondences (in the procedures
Associate() andResetLabelMap() )suchthat
each entry is assignedthe value of the associatedfinal
label. Subsequently, the labels of the runs are updatedusing this LabelMap . At last, the labeled image can
be written as a pixel array by decoding the run length
representation RleData .
In the pseudo-code for the labeling method
DoLabeling() , variables which are not described
belowhavethesamemeaningasin DoEncoding() .
NB_VALID – Boolean variable which indicates
whether a neighbor pixel lies inside the image
(true)ornot( false)
newPrelLabel – variable for assigning a unique
preliminary label to each non-labeled run (with
placeholderlabel)if LABEL_PROP == true
prelLabel – stores the preliminary label of a run;
ifLABEL_PROP == false , thenprelLabel
has the same meaning as in DoEncoding() ,
otherwise prelLabel does not necessarily
increase while passing through the image but
mightholdapropagatedlabel
CurNbCs – vector holding the coordinates of the
current neighborpixel (w.r.t. the original order of
directions)
LabelMap [global]–vectormappingthepreliminary
labels (as the indices of the entries) to the
associatedfinallabels(astheirvalues)
The labeling procedure DoLabeling() returns
the number of objects in Imagecorresponding to
the adjacency deducible from the predefined array
Neighbors . In the following pseudo-code, code
passages enclosed by ‘ ###’ should be regarded as
present only if the condition in the ###-header is
fulfilled.
int DoLabeling()
{
bool NB_VALID;
int i, j, k, m, r, newPrelLabel, prelLabel, objects;
vector<int> CurPixCs(n), CurNbCs(n), LinePos(n-1), Labe lMap;
RleLine *pLine;
54
ImageAnalStereol2009;28:45-61
RleRun *pRun, *pNbRun;
// set variables
newPrelLabel = 0;
CurPixCs.init(0);
LinePos.init(0);
IF LABEL_PROP THEN
LabelMap.resize(1);
ELSE
LabelMap.resize(prelLabel+1);
LabelMap.init(BG_COL);
// main loop
DO
{
pLine = &RleData[LinePos[n-2]]..[LinePos[0]];
// set the coordinates of the current pixel
FOR r FROM 1 TO n-1 DO
CurPixCs[Index[r]] = LinePos[r-1];
// test all runs in the current line addressed by pLine
FOR i FROM 0 TO pLine->size()-1 DO
{
pRun = &(pLine->at(i));
### LABEL_PROP == true
// check whether run has yet received a propagated label
// otherwise assign a new preliminary label and register it
IF pRun->label == -999 THEN
{
pRun->label = ++newPrelLabel;
LabelMap.add(BG_COL);
}
###
prelLabel = pRun->label;
// examine the start and endpoint of the current run
FOR j FROM 0 TO 1 DO
{
IF j == 0 OR (j == 1 AND pRun->Pos[0] != pRun->Pos[1]) THEN
{
CurPixCs[Index[0]] = pRun->Pos[j];
// test the valid neighbors for label correspondences
FOR k FROM 0 TO Neighbors.size()-1 DO
{
NB_VALID = true;
FOR m FROM 0 TO n-1 WHILE NB_VALID DO
{
CurNbCs[m] = CurPixCs[m] + Neighbors[k][m];
IF CurNbCs[m] < 0 OR CurNbCs[m] >= ImgSize[m] THEN
NB_VALID = false;
}
// neighbor is valid, i.e. lies inside the image
IF NB_VALID THEN
{
55
SANDFORT KET AL:Labelingofn-dimensionalimages
// query the address of the current neighbor
// this is NULL if the neighbor is a background pixel
pNbRun = QueryPtr(CurNbCs);
IF pNbRun != NULL AND pNbRun->label != prelLabel THEN
{
### LABEL_PROP == true
IF pNbRun->label != -999 THEN
Associate(prelLabel, pNbRun->label, LabelMap);
ELSE
pNbRun->label = prelLabel;
###
### LABEL_PROP == false
Associate(prelLabel, pNbRun->label, LabelMap);
###
}
}
}
}
}
}
// update the line position, according to the scanning order
r = 0;
WHILE r < n-1 AND ++LinePos[r] >= ImgSize[Index[r+1]] DO
LinePos[r++] = 0;
}
WHILE LinePos[n-2] < ImgSize[Index[n-1]]
// setup the final LabelMap
objects = ResetLabelMap(LabelMap);
// update the run labels
LinePos.init(0);
DO
{
pLine = &RleData[LinePos[n-2]]..[LinePos[0]];
FOR i FROM 0 TO pLine->size()-1 DO
pLine->at(i).label = LabelMap.at(pLine->at(i).label);
r = 0;
WHILE r < n-1 AND ++LinePos[r] >= ImgSize[Index[r+1]] DO
LinePos[r++] = 0;
}
WHILE LinePos[n-2] < ImgSize[Index[n-1]]
// return the number of objects in the image
return objects;
}
The array RleData now contains all information
on the labeled image. The next routine decompressesthis information andwrites animagein pixelformat.
void WriteLabeledImage()
{
int i, r;
vector<int> CurPixCs(n), LinePos(n-1);
RleLine *pLine;
RleRun *pRun;
56
ImageAnalStereol2009;28:45-61
LinePos.init(0);
DO
{
pLine = &RleData[LinePos[n-2]]..[LinePos[0]];
FOR r FROM 1 TO n-1 DO
CurPixCs[Index[r]] = LinePos[r-1];
FOR i FROM 0 TO pLine->size()-1 DO
{
pRun = &(pLine->at(i));
FOR CurPixCs[Index[0]] FROM pRun->Pos[0] TO pRun->Pos[1] DO
Image[CurPixCs[0]]..[CurPixCs[n-1]] = pRun->label;
}
r = 0;
WHILE r < n-1 AND ++LinePos[r] >= ImgSize[Index[r+1]] DO
LinePos[r++] = 0;
}
WHILE LinePos[n-2] < ImgSize[Index[n-1]]
}
The overall calling sequence is DoEncoding()
-DoLabeling() -WriteLabeledImage() .
We now explain shortly the auxiliary routines
QueryPtr() ,Associate() , as well as
ResetLabelMap() .
Themethod QueryPtr() returnsapointertothe
runwhichcontainsthepixelwithcoordinatevector Cs
if this is not a background pixel. Otherwise, the pixel
is not contained in any run, and the NULL-pointer is
returned.
RleRun *QueryPtr(vector<int>& Cs)
{
int i, pos;
vector<int> LinePos(n-1);
RleLine *pLine;
pos = Cs[Index[0]];
FOR i FROM 1 TO n-1 DO
LinePos[i-1] = Cs[Index[i]];
pLine =
&RleData[LinePos[n-2]]..[LinePos[0]];
FOR i FROM 0 TO pLine->size()-1 DO
{
IF pLine->at(i).Pos[1] >= pos THEN
{
IF pLine->at(i).Pos[0] <= pos THEN
return &(pLine->at(i));
ELSE
return NULL;
}
}
return NULL;
}The recursive method Associate() relates
equivalent preliminary labels aandbwith the
smallest equivalent label in LabelMap . After calling
this routine, traversing LabelMap via repeated
substitution of indexbyLabelMap[index] as
long as the latter is non-zero, starting from either
index = a orindex = b , leads to the smallest
equivalentlabelasthe last index.
void Associate(int a, int b,
vector<int>& LabelMap)
{
int c;
IF a < b THEN
{ c = a; a = b; b = c; }
IF LabelMap[a] == BG_COL THEN
LabelMap[a] = b;
ELSE IF LabelMap[a] != b THEN
Associate(LabelMap[a], b, LabelMap);
}
Finally, ResetLabelMap() assigns unique
and minimal final labels to the preliminary ones
which have no smaller counterparts and resolves the
connectionsmade up by Associate() .This means
that afterwards each preliminary label iis mapped
to the smallest integer j = LabelMap[i] which
preserves the equivalence relations. The return value
ofResetLabelMap() is the number of objects in
Image.
int ResetLabelMap(vector<int>& LabelMap)
{
int i, objects;
57
SANDFORT KET AL:Labelingofn-dimensionalimages
objects = 0;
FOR i FROM 0 TO LabelMap.size()-1 DO
{
IF LabelMap[i] == BG_COL THEN
LabelMap[i] = ++objects;
ELSE
LabelMap[i] = LabelMap[LabelMap[i]];
}
return objects;
}
The latter both routines have been proposed
and discussed in Mart´ ın-Herrero and Pe´ on-Fern´ andez
(2000). The algorithm presented in this section
is the basis for the current labeling function in
the commercial software MAVI (Fraunhofer ITWM,
DepartmentofImageProcessing,2005).
EXAMPLEANDDISCUSSION
In this final section, we look at an exemplary
application of our algorithm as well as of the
alternative recursive algorithm by Mart´ ın-Herrero,
cf.Mart´ ın-Herrero (2004). The test data are two
binarized 3D microtomography images of materials.
Theirvisualizationsareshownin Figs.1 and2.
Fig.1.3Dimage(660 ×660×660pixels,0.7µmpixel
size) showing a specimen of a metallic foam in early
extensionstage.
InFig.1themetallicmatrixistransparentwhilethe
porespaceappearsopaque.InFig.2thesolidmatteris
opaqueandthe pore spaceis transparent.Bothimagesare at hand in RAW format with 8 Bit per pixel. The
firstimagehasasizeof274MBand5.7%foreground
pixels. The second one needs 26 MB and has 86.4%
foreground pixels. For details about the materials and
imagingmethodsseeHelfen etal.(2002;2003).
Our test system is a PC with AMD Athlon(tm) 64
Processor 3800+, 1 GB RAM, running SuSE Linux
10.0 (kernel v2.6.13). We compiled the source codes
with the GCC v4.0.2, using the optimization switch
‘-O3’. The performance data are listed in the tables
below (‘S./O.’ denotes our algorithm, ‘M.-H.’ the one
from Mart´ ın-Herrero).
Fig.2.3Dimage(299 ×300×300pixels,0.75µmpixel
size) showing a specimen of the Fontainebleau sand
stone.
At first, we want to make clear what the
times listed in Tables 1 and 2 refer to. For our
algorithm, we note that after calling the procedures
DoEncoding() andDoLabeling() the run
length array RleData with the correct run labels is
obtained.Itmightbesufficientorevenbeneficialtodo
any further processing of the image on this structure.
Otherwise, it remains to write the labeled image as
a pixel array using WriteLabeledImage() . The
processing time for this has not been considered in
Tables1and2.SincethealgorithmbyMart´ ın-Herrero
identifies non-labeled foreground pixels by a negative
value, it usually necessitates a preprocessing of the
source pixel array Imageto change the foreground
labelFG_COL to this value (if the values in Image
arenon-negative).In principle,theidentificationvalue
can be any value which does not coincide with the
58
ImageAnalStereol2009;28:45-61
background color BG_COL (mostly chosen to be 0)
and any valid final label. However, the number of
objects in Imageis not predicted, and we require the
final labels to be consecutive integers starting from 1.
If one agrees to drop this restriction, a substitution of
FG_COL can be avoided. Anyway, the time possibly
needed for this is not recorded above. We also want
to emphasize that the times have been obtained with
variants of both algorithms optimized for 3D data and
a fixedscanningorder.
In comparison with the second image, the first
one has only few foreground pixels and a pixels per
object ratio of 59. By contrast, the second image
containsmainlyforegroundpixels,whichare multiply
connected, and has 704926 pixels per object in
average. On the first dataset, the recursive algorithm
by Mart´ ın-Herrero is more efficient than ours and
takes clear advantage of extremal label propagation.
Nevertheless,itappearsthatourmethodwithactivated
label propagation can serve as an alternative. Since
the algorithms rely on different image representations
(pixel array vs.RLE structure), it is not appropriate
to compare the times for the labeling only. In this
respect, we remark that it could be valuable for both
algorithms to inspect how they can be adapted to
different image representations. For this question, the
fundamental work Dillencourt et al.(1992) should be
consulted.
On the dataset for Fig.2, our approach combined
with label propagation shows some superiority
compared to the recursive one. For the latter, we
note here the relatively high recursion depth along
with the fact that many pixels are queried which
have already received a propagated label. The label
propagation in our method significantly reduces the
numberofpreliminarylabelsforthisimageandmakes
the labeling procedure including the computation
ofLabelMap very fast. For the evaluation andcomparison, it is also important that the runs for the
second image contain 39 pixels in average, while for
the first one this number is 4. Roughly speaking,
it appears that a setting in which many pixels are
captured in few runs overlapping in adjacent lines is
more amenable to our method than to the other. Also
in view of the above times, it should be carefully
observed whether a preprocessing such as RLE leads
to an overall improvement in the processing needs
(the processing time, in particular). Although for
the worst case of a ‘salt and pepper’ image (with
a large amount of small objects) also the recursive
and probably almost every other labeling algorithm
can not fully deploy its potential, a run length
encodingsurelywill spoilthe totalperformanceofthe
procedure. Mart´ ın-Herrero’s algorithm from Mart´ ın-
Herrero (2004) and the enhanced Hoshen-Kopelman
algorithm from Hoshen (1998) offer an on-the-fly
analysis of objects. This feature could be added to
our algorithm by collecting information during the
run length encoding and processing it subject to the
label correspondences (and propagations) during the
labeling.
Twootherimportantimpactfactorsforthelabeling
are the adjacency system and the scanning order.
The above results concern the labeling w.r.t. the 6-
adjacency. Since neighbors in the direction of the
(quasi-)run length encoding need not to be checked
for correspondences, a κ-adjacency requires a test of
onlyκ−2 neighbors in DoLabeling() and the
recursion part of the algorithm from Mart´ ın-Herrero
(2004), respectively.When switching to e.g.the 14.x-
adjacency, in our method the ratio (14−2)/(6−
2) =3 of tested neighbors is roughly reflected in
the times needed for the labeling. Although described
onlyforthe6-adjacencyinMart´ ın-Herrero(2007),the
recursivemethodcanbestraightforwardlyextendedto
otheradjacencies.A testandcomparisonforthesehas
notyetbeendone.
Table1.Performancedatafor Fig.1andthe6-adjacency,275536obje cts.
algorithmS./O.: ♯preliminary labels
M.-H.: max. recursiondepthexecutiontime
RLELabeling total
S./O.,LABEL_PROP = false 3856027 2.73s5.93s8.66s
S./O.,LABEL_PROP = true 738601 2.72s3.70s6.42s
M.-H. 504706 —4.37s4.37s
Table2.Performancedatafor Fig.2andthe6-adjacency,33 objects.
algorithmS./O.: ♯preliminary labels
M.-H.: max. recursiondepthexecutiontime
RLELabeling total
S./O.,LABEL_PROP = false 590753 0.34s2.96s3.3s
S./O.,LABEL_PROP = true 28176 0.33s0.58s0.91s
M.-H. 590649 —1.83s1.83s
59
SANDFORT KET AL:Labelingofn-dimensionalimages
Regarding the theoretical foundation from section
“Connectedness”, we want to remark the following.
Both algorithms clearly work for arbitrary adjacency
systems including the 18-adjacency. However, one
should be aware that further processing and analysis
of images labeled w.r.t. the 18-adjacency can
lead to conflicts. Since there does not exist an
adjacency system for the background which is
complementary to the 18-adjacency according to
Definition 3, the topology of the foreground pixels
does not fit the topology of the background pixels.
In particular, the Euler number of the background
differs from the sum of the Euler numbers of the
obtained equivalence classes. Thus, whenever further
processingandanalysisofthelabelimageisnecessary,
we recommend a labeling w.r.t. an adjacency system
Fwith existing Fcaccordingto Definition 3.
A change in the scanning order can have
significant impact on the processing time, depending
on the anisotropy of the image. Using the general
implementationsofthealgorithms(notoptimizedfora
specificscanningorder),thechangeinprocessingtime
for the 6-adjacency is as follows: For the algorithm
from Mart´ ın-Herrero (2004), it is up to 11% for the
first image and up to 13% for the second one. For
our algorithm with label propagation, it is up to 19%
for the first image and up to 7% for the second
one. However, one should regard the short times
and a possible amplification of marginal variations.
Moreover, the second image is nearly isotropic. The
impact of the choice of the scanningorder is expected
to become obvious for strongly anisotropic images.
For the datasets used in this benchmark and further
detailsseethewebpageathttp://www.itwm.fhg.de/bv/
projects/BA/RLE.
CONCLUSION
In the first part of our paper, we discussed the
discretization and labeling of an n-dimensional set
w.r.t. to a given adjacency system. We suggest here
to define the complementarity of adjacency systems
by a relation for the Euler number instead of by the
Jordan-VeblenandJordan-Brouwertheoremin2Dand
3D, respectively.This appearsnatural from a practical
pointofview andis relevantforthe consistencyofthe
labelingresults.Inthesecondpart,weproposedanew
and flexible labeling algorithm based on a preceding
run length encoding of the input image. We shortly
addressheretheapplicationsofthispreprocessingand
a fewpossibleextensionsofthe algorithm.
If the image or its objects are to be further
processedafter the labeling, then a compactand easy-
manageable representation of the image as given in asimpleformbyarunlengtharraycanbedesirable.We
only mention the morphological filtering of objects,
noise removal and geometric transforms as operations
which can take advantage of this. For further uses
of the run length encoding we refer to the papers
Di Zenzo et al.(1996) and Messom et al.(2002).
SincetherecursivemethodofMart´ ın-Herreroandours
follow in a way complementary ideas, a combination
ofbothmethodscouldbeworth considering.
Moreover, the tasks of run length encoding and
labeling are well-suited for parallelization. A simple
step towards this would be to separatebig images into
blocks and process first their interiors independently
and then gather correspondences in the border areas.
Using label offsets after the block computations to
avoida falsecoincidenceoflabels,a labelmap forthe
wholeimagecanbecomputed.
ACKNOWLEDGMENTS
The authors very much appreciate the comments
and suggestions by the referees. We also thank Bj¨ orn
Wagner and Michael Godehardt from the Fraunhofer
Institute for Industrial Mathematics (ITWM) for
software support and fruitful discussions about
implementation of algorithms. The research of the
second author was supported by the FH3-programme
of the German Federal Ministry of Education and
Researchunderprojectgrant1711B06.
REFERENCES
Aguilera A, Rodr´ ıguez J, Ayala D (2002). Fast connected
component labeling algorithm: A non voxel-based
approach. Tech. Rep. LSI-02-26-R, Universitat
Polit` ecnica, Catalunya. http://www.lsi.upc.edu/dept/
techreps/.
ArmstrongMA (1997). Basic Topology. Berlin: Springer.
Di Zenzo S, Cinque L, Levialdi S (1996). Run-based
algorithms for binary image analysis and processing.
IEEETransPatternAnal18:83–9.
Dillencourt MB, Samet H, TamminenM (1992). A general
approachtoconnected-componentlabelingforarbitrary
imagerepresentations. JACM 39:253–80.
Fraunhofer ITWM, Department of Image Processing
(2005). MAVI – Modular Algorithms for Volume
Images. http://www.itwm.fhg.de/mab/projects/MAVI/.
Helfen L, Baumbach T, Schladitz K, Ohser J (2003).
Determinationofstructuralpropertiesoflightmaterials
by three-dimensional synchrotron-radiation imaging.
ImagingMicrosc5:55–7.
Helfen L, Baumbach T, Stanzick H, Banhart J,
Elmoutaouakkil A, Cloetens P (2002). Viewing
the early stage of metal foam formation by computed
60
ImageAnalStereol2009;28:45-61
tomography using synchrotron radiation. Adv Eng
Mater4:808–13.
HopcroftJ,TarjanRE(1973).Efficientalgorithmsforgraph
manipulation. CommACM16(6):372–8.
Hoshen J (1998). On the application of the enhanced
Hoshen-Kopelman algorithm for image analysis.
PatternRecognLett 19:575–84.
Hoshen J, Koppelman R (1976). Percolation and cluster
distribution. I. cluster multiple labeling technique and
critical concentration algorithm. Phys Rev B 14:575–
84.
Klette R, Rosenfeld A (2004). Digital Geometry.
Amsterdam:Morgan&KaufmanPubl.
Lachaud JO, Montanvert A (2000). Continuous analogs
of digital boundaries: A topological approach to iso-
surfaces. GraphicalModels62:129–64.
Mart´ ın-Herrero M (2004). Hybrid cluster identification. J
PhysA MathGen37:9377–86.
Mart´ ın-HerreroM(2007). Hybridobjectlabellingindigit al
images. MachineVisionAppl18:1–15.
Mart´ ın-Herrero M, Pe´ on-Fern´ andez J (2000). Alternativ e
techniquesforcluster labellingon percolationtheory. J
PhysA MathGen33:1827–40.
Messom CH, Demidenko S, Subramaniam K, Gupta GS
(2002). Size/position identification in real-time image
processing using run length encoding. Proc 19th IEEE
InstrumMeasTechConf2:1055–9vol.2.
Nagel W, Ohser J, Pischang K (2000). An integral-
geometric approach for the Euler-Poincar´ e
characteristicofspatial images. JMicrosc198:54–62.
Ohser J, Nagel W, Schladitz K (2002). The Euler number
of discretized sets – on the choice of adjacency inhomogeneous lattices. In: Mecke KR, Stoyan D, eds.,
Morphology of Condensed Matter, Lect Notes Phys.
Berlin:Springer.
OhserJ,NagelW,SchladitzK(2003). TheEulernumberof
discretisedsets–surprisingresultsinthreedimensions.
ImageAnalStereol22:11–9.
Ritter GX, Wilson JN (2001). Handbook of computer
vision algorithmsin image algebra,chap. 6, Connected
Component Algorithms. Boca Raton, Florida: CRC
Press, 173–86.
RosenfeldA(1970). Digitaltopology. AmerMathMonthly
86:621–30.
Rosenfeld A, Pfaltz JL (1966). Sequential operations in
digitalpictureprocessing. JACM13:471–94.
Rotman JJ (1993). An Introductionto Algebraic Topology.
Berlin:Springer-Verlag.
Schladitz K, Ohser J, Nagel W (2006). Measurement of
intrinsic volumes of sets observed on lattices. In:
Kuba A, Nyul LG, Palagyi K, eds., Proc 13th Int
Conf Discrete Geom Comput Imagery, LNCS. DGCI,
Szeged,Hungary,Berlin:Springer.
Schmitt M (1998). Digitization and connectivity. In:
Heijmans HJAM, Roerdnik JBTM, eds., Mathematical
morphology and its application to image and signal
processing.Dortrecht:KluwerAcademicPublishers.
SchneiderR(1993).ConvexBodies:TheBrunn-Minkowski
Theory. Cambridge: Encyclopediaof Mathematicsand
ItsApplicationVol. 44,CambridgeUniversityPress.
Thurfjell L, Bengtsson E, Nordin B (1992). A new three-
dimensional connected components labeling algorithm
with simultaneous object feature extraction capability.
ComputVisionGraph54:357–64.
61