Image Anal Stereol 2016;35:167-179 doi:10.5566/ias.1489
Original Research Paper
ESTIMATINGFIBREDIRECTIONDISTRIBUTIONSOFREINFORCED
COMPOSITESFROMTOMOGRAPHICIMAGES
OLIVER WIRJADI
1
, KATJA SCHLADITZ
B
;1
, PRAKASH EASWARAN
1;2
AND JOACHIM
OHSER
3
1
Fraunhofer ITWM, Fraunhofer-Platz 1, D-67663 Kaiserslautern, Germany;
2
University of Kaiserslautern,
Germany;
3
Hochschule Darmstadt, Fachbereich Mathematik und Naturwissenschaften, Sch¨ offerstraße 3,
D-64295 Darmstadt, Germany
e-mail: oliver.wirjadi@gmail.com, katja.schladitz@itwm.fraunhofer.de,
prakash.easwaran@itwm.fraunhofer.de, joachim.ohser@h-da.de
(Received February 24, 2016; revised October 10, 2016; accepted November 2, 2016)
ABSTRACT
Fibre reinforced composites constitute a relevant class of materials used chiefly in lightweight constructions
for example in fuselages or car bodies. The spatial arrangement of the fibres and in particular their direction
distribution have huge impact on macroscopic properties and, thus, its determination is an important topic
of material characterisation. The fibre direction distribution is defined on the unit sphere, and it is therefore
preferable to work with fully three-dimensional images of the microstructure as obtained, e.g., by computed
micro-tomography. A number of recent image analysis algorithms exploit local grey value variations to
estimate a preferred direction in each fibre point. Averaging these local results leads estimates of the
volume-weighted fibre direction distribution. We show how the thus derived fibre direction distribution is
related to quantities commonly used in engineering applications. Furthermore, we discuss four algorithms
for local orientation analysis, namely those based on the response of anisotropic Gaussian filters, moments
and axes of inertia derived from directed distance transforms, the structure tensor, or the Hessian matrix.
Finally, the feasibility of these algorithms is demonstrated for application examples and some advantages and
disadvantages of the underlying methods are pointed out.
Keywords: composites, fibre direction distribution, image analysis, orientation tensor, tomography.
INTRODUCTION
Light-weight construction in various contexts
is one of the main driving forces behind the
increasing application of fibre-reinforced materials.
The microstructure of these materials is largely
responsible for their favorable mechanical properties.
Especially the directional distribution of the fibre
system influences their tensile strength. This fact is
also accounted for in modern tools for macroscopic
material property simulation, leading to an increasing
demand for precise measurement of local fibre volume
content and fibre directions.
Early work has focused on the estimation of fibre
direction distributions from planar sections through
the material. Obviously, a planar cross-section of a
straight cylinder is up to edge effects an ellipse whose
aspect ratio depends on the angle between section
plane and cylinder core. This observation was used to
estimate the fibre direction, e.g., by Zhu et al. (1997).
However, a fibre direction is not uniquely determined
by the corresponding cross-section. To overcome this
problem, Clarke et al. (1995) and Eberhardt and Clarke
(2001) used confocal laser scanning microscopy to
determine fibre directions by tracking individual fibres
in glass fibre reinforced polymers through a depth of
up to 150 mm. Nevertheless, this type of stereological
reconstruction of the three-dimensional fibre directions
from planar sections is error-prone in particular for
curved fibres or for straight fibres which are parallel
to the section plane. These problems can be avoided
by computed micro-tomography (mCT), which is
capable to fully reconstruct the three-dimensional
microstructure of fibre reinforced materials.
Recently, a number of quantities characterising
fibre direction distributions as well as a wide range
of image analysis algorithms for computing these
quantities from 3D data obtained by mCT have been
proposed. The rose of directions, i.e., the probability
density function of the direction distribution, can be
obtained via the inverse cosine transform of data
either obtained from so-called generalised projections
onto subspaces (Ohser and Schladitz, 2009) or from
counts of fibre sections in section planes (Kiderlen and
Pfrang, 2005; Riplinger and Spiess, 2012). While the
former method is highly efficient with a low directional
and high lateral resolution, the latter one features a
high directional resolution and a low lateral resolution.
Ohser et al. (1998) pointed out the mutual exclusion
of increasing both, directional and lateral resolution
167

WIRJADI O ET AL: The volume-weighted direction distribution
for this class of algorithms. In the present paper,
estimating the rose of directions is not considered,
since in contrast to the methods that we will use here,
it does not allow for a local analysis of fibre systems.
The segmentation of fibre profiles either in section
planes is a non-trivial image processing task. In order
to compute fibre direction distributions without this
tedious pre-processing, some authors have proposed to
reduce the computation of fibre direction distributions
to the computation of local fibre directions in
individual pixels. As all of these methods average over
pixels to obtain local directional estimates, they can
be characterised as algorithms for a volume-weighted
fibre direction estimation.
Robb et al. (2007), Wirjadi (2009), and Wirjadi
et al. (2009) test for each fibre pixel the filter response
of anisotropic, prolate Gaussian filters with major
principal axis from a finite set of directions, see the
section “Methods” below for details. As we will show,
such an approach is viable but of limited accuracy
due to the strongly restricted number of directions
that can be considered in practice due to runtime
limitations. Sandau and Ohser (2007) proposed the so-
called chord length transform (CLT) which assigns to
each fibre point the chord of maximum length covered
by the fibre. The direction of this chord is used as
an estimate of the local fibre direction. The CLT will
not be considered within the present paper. It requires
testing a certain number of chord directions, entailing
the same problematic dependence between accuracy
and runtime as observed for the anisotropic Gaussian
filter method. The approach of Altendorf and Jeulin
(2009) can be seen as an advanced CLT as it is based
on the CLT idea but refines estimation of the local fibre
direction by computing it from the main axis of inertia.
This latter approach will be summarised and evaluated
below.
Choosing a finite number of test directions is not
required at all when considering partial derivatives of
the local grey value structure along fibres. Krause et al.
(2010) accumulate gradient information in structure
tensors and derive the local fibre direction from the
eigendecomposition of these tensors. Similarly, the
Hessian matrix of second partial derivatives yields
the local fibre direction, as this corresponds to the
direction of lowest grey value curvature as observed
by Frangi et al. (1998); Tankyevych et al. (2009). Both
of these approaches based on partial derivatives will
be described in more detail in the following section
“Methods”.
This paper contributes to the development of
fibre direction estimation by subsuming the previously
proposed algorithms in a common framework and
by showing how the volume-weighted fibre direction
distribution relates to other quantities characterising
directional features of fibre structures. Lastly, our
numerical results on simulated data may be used as
a guideline for practitioners seeking suitable analysis
algorithms.
MATERIALSANDMETHODS
ORIENTATIONANALYSIS
In this section, we introduce methods that allow
us to characterise the directional distributions of
fibre systems in composite materials based on image
data. In particular, we define the volume-weighted
directional distribution, derive the orientation tensors,
and describe the four methods for estimating the fibre
direction in each pixel. LetF denote a macroscopically
homogeneous (but anisotropic) random fibre system
in R
3
following the definition by Mecke and Nagel
(1980); Nagel (1983), that isF is a random system of
rectifiable curves.
In our setting, the thickness of the fibres is not
negligible. We will therefore address the random set
X =F B
r
as the fibre system while calling F the
fibre core system. That is, the random fibre system
X is derived from the random fibre core system F
by dilation with the ball B
r
of radius r > 0. Note
that for a wide range of fibre-reinforced composites,
e.g., most glass and carbon fibre-reinforced polymers
and steel fibre-reinforced concretes, the fibre radius
r is known from the production process. Therefore,
we will assume fibres with a circular cross-section
and constant diameter 2r throughout this paper. The
fibre system X may be sufficiently smooth, and
self-intersections may not occur. More precisely, we
suppose that there is an e > 0 such that X is
morphologically closed w. r. t. B
e
,
(X B
e
) B
e
=X: (1)
Condition Eq. (1) ensures that the direction of the fibre
core system is uniquely defined in each point within
the fibre system. Moreover, being morphologically
closed means that the local curvature of the fibre core
system must be less than 1=r and thus in particular that
the fibre surfaces are not rough.
Finally, assume that fibres are fibres in the sense
that they are longer than wide, that is, their length is
significantly larger than their diameter 2r.
For a point y in the fibre core system F, denote
by v(y) the tangent direction. We now define the local
fibre direction v(x) of a point x in the fibre systemX.
168

Image Anal Stereol 2016;35:167-179
This direction corresponds to the fibre core system’s
tangent direction in the point that is closest to x, i. e.,
v(x)= v(y) with y= argmin
y2F
kx yk. This point y
is unique due to the assumption Eq. (1) onF and the
diameter of the fibres in X made above. The vectors
v(x) form a random vector field. Note that we need not
to consider the orientation, i. e., the sign of v(x), for
its choice in a fibre point x2X is ambiguous. Thus we
consider the set S
2
+
of non-oriented directions that is
equivalent to the upper hemisphere.
Disregarding noise and artefacts for the moment,
a tomographic reconstruction yields a 3D image f
which is an observation of the fibre system X of the
form f(x) = 1
X
(x), where 1
X
denotes the indicator
function of X. Now, assume for the moment that we
were able to compute the random vector field v(x) at
every point x2X in the fibre system from an image
f . In the current context, we are interested in the fibre
direction distribution in the typical fibre point of the
fibre system. Then the mean R over local directions
v(x) in a compact window W  R
3
with volume
vol(W)> 0,
R(A)=
1
2pEvol(X\W)
E
Z
X\W
1
A
(v(x))dx; (2)
where A S
2
+
is a measurable set of directions, is a
probability measure. The expectation in Eq. (2) should
be understood with respect to the fibre systemX. The
distribution R is volume-weighted, since it averages
over all points occupied by the fibre systemX.
For systems of fibres with variable thickness,
the just defined direction distribution R is a so-
called volume-weighted one. More precisely, R(A)
is the direction distribution function of the typical
point ofX. Analogously, a surface-weighted direction
distribution is the direction distribution of the typical
point on the surface¶X ofX and, finally, the direction
distribution of the core system F is the length-
weighted direction distribution as defined by Mecke
and Nagel (1980). These distributions can differ
considerably. In our setting of constant fibre thickness
however, the volume-, surface-, and length-weighted
fibre direction distributions are identical (up to edge
effects resulting from contributions of the fibre ends).
To summarise, if we were given access to the
random vector field v(x) and to the set of points
occupied by the fibre system X in an observed and
possibly degenerated image f , the volume-weighted
fibre direction distribution would be available with
Eq. (2).
ORIENTATIONTENSORS
Instead of the probability measure R, it is often
more convenient to work with the moments of R, the
so-called orientation tensors as defined by Advani and
Tucker (1987). Let u
i
denote the ith component of
some normalsed direction vector u. Then the second
and fourth order orientation tensors are defined as
A
(2)
=(a
i j
) and A
(4)
=
 a
i jk‘
 with
a
i j
=
Z
S
2
+
u
i
u
j
R(du);
a
i jk‘
=
Z
S
2
+
u
i
u
j
u
k
u
‘
R(du); i; j;k;‘= 1;2;3:
Tensors of orders other than two and four can of course
be defined in a similar way. In order to highlight
the relevance of orientation tensors, we follow the
arguments of Advani and Tucker (1987): Material
properties such as the behaviour under mechanical
stress are usually modeled in continuum, where a
continuous, homogeneous material is assumed, rather
than a two-phase medium such as fibre-reinforced
polymers. This process requires averaging of the
properties of the two materials over the direction
distribution R. The average of a material property T
over all space directions is denoted byhTi,
hTi=
Z
T(u)R(du): (3)
For computational efficiency, it is desirable to replace
the integral over the hemisphere in Eq. (3) by some
simplified expression. This is where orientation tensors
fit in. For instance, under some symmetry assumptions
for the material (called transverse isotropy by Advani
and Tucker (1987)), the second-order, tensor-valued
material property tensor T
(2)
is a sufficient description
ofhTi, and can be decomposed as
T
(2)
=(T
i j
) with T
i j
= b
1
a
i j
+ b
2
d
i j
;
where b
1
and b
2
are scalars, and d
i j
is the Kronecker
delta. Depending on the symmetries, similar relations
can also be derived between fourth-order tensor-valued
properties and the fourth-order orientation tensor A
(4)
.
Examples for second-order tensor-valued properties
are thermal expansion and thermal conductivity, while
elastic stiffness or viscosity can be described as fourth-
order tensors (Advani and Tucker, 1987; Bunge, 1993).
Thus, orientation tensors play a central role in the
simulation of material properties.
Beyond their use in continuum modeling,
orientation tensors, especially the eigendecompo-
sition of A
(2)
, allow for intuitive compact descriptions
of the shape of the direction distribution R. For
169

WIRJADI O ET AL: The volume-weighted direction distribution
instance, the eigenvector to the largest eigenvalue
l
max
of A
(2)
corresponds to the mean fibre direction.
Furthermore, ratios between a
i j
’s eigenvalues can be
used to characterise the anisotropy of R, see Fisher
et al. (1987). One such useful characteristic is the
anisotropy factor a = 1 l
min
=l
max
. This degree of
anisotropy a is zero for the isotropic fibre direction
distribution. a = 1 means that the fibre direction
distribution is concentrated on a great circle. In the
latter case, a does however not differentiate between
a uniaxial or an isotropic distribution within the plane
determined by the great circle.
These two scenarios motivate our search for
methods that compute the orientation tensors A
(2)
and
A
(4)
directly from an image f . We follow the same
concept as for the directional fibre distribution in
Eq. (2) and replace the integral w. r. t. the unknown
probability measure R in Eq. (3) by a sample over local
fibre direction vectors v(x) with components v
i
(x).
This yields
a
i j
=
1
Evol(X\W)
E
Z
X\W
v
i
(x)v
j
(x)dx; (4)
a
i jk‘
=
1
Evol(X\W)
E
Z
X\W
v
i
(x)v
j
(x)v
k
(x)v
‘
(x)dx:
(5)
for the components of the second and fourth order
orientation tensors.
The remainder of this paper is dedicated to the
introduction of algorithms for estimating R and the
derived characteristics from image data of fibre-
reinforced composites, to their interpretations, and to
the experimental evaluation and application of those
methods.
IMAGEBINARISATION
From an image processing perspective, the
quantities introduced above simplify the problem at
hand. Methods for computing the fibre direction
distribution proposed so far usually required precise
segmentations of fibres either in 2D (Clarke et al.,
1995; Lee et al., 2002; Kiderlen and Pfrang, 2005),
or in 3D (Yang and Lindquist, 2000; Tan et al.,
2006). In contrast to that, Eq. (2) and Eq. (5) merely
require the fibre system X. That is, these estimators
do not require a segmentation of individual fibres, but
only a segmentation of the complete fibre system. A
binarisation of an image is a separation of its pixels
into two disjoint sets. In the present context, it yields
an estimation of the set of pixels covered by the fibre
system and its complement.
This can be achieved using standard image
processing methods. From our experience, sufficient
binarisations of tomographic reconstructions of fibre
reinforced polymers can easily be achieved using some
low-pass noise reduction filters and a global grey value
threshold. To demonstrate that this is true for images
encountered in practice, consider the binarisation of
the mCT image of a glass fibre-reinforced polymer
in Fig. 1. Here, a simple smoothing filter (3 3 3
discrete Gaussian) and a global threshold value derived
from the grey value distribution in the image were
used. While such a binarisation is not sufficient to
segment individual fibres, it is sufficient as an estimate
of the fibre systemX required in the present paper.
(a) Original (b) Binarisation
Fig. 1. Binarisation of the fibre system of a glass fibre-
reinforced polymer (GFRP) imaged at 3mm pixel size.
Further results regarding this dataset will be presented
in the Application section below. Here, a 3 3 3
discrete Gaussian filter is followed by a simple global
threshold.
Note that the thus derived binarisation might
locally not be good enough to fulfill the condition
Eq. (1). Thus, some pixels of the fibre system can not
be assigned uniquely to a fibre core. However, as this
problem typically occurs where parallel fibres touch, it
does not cause large estimation errors.
For all algorithms introduced below, we will use
such image binarisations as an estimation ofX in order
to restrict the evaluation of fibre directions to pixels
belonging to the fibre system.
ESTIMATIONOFLOCALFIBRE
DIRECTIONS
Local directional information has widely been
used in the image processing literature. The moments
of inertia have been used as features for object and
pattern recognition for a very long time, see e.g., Hu
(1962); Prokop and Reeves (1992), and more recently
to analyse and segment images of fibrous materials by
Altendorf and Jeulin (2009; 2011) and Altendorf et al.
(2012).
170

Image Anal Stereol 2016;35:167-179
Anisotropic linear filters, and the directional
information contained in first and second order
derivatives have been exploited for image filtering:
Anisotropic Gaussian filters were shown to be suitable
for adaptive smoothing of diffusion tensor magnetic
resonance data (DT-MRI), e.g., by Lee et al. (2006).
The structure tensor (a 3 3 matrix of smoothed
pairwise products of first partial derivatives) is used
in diffusion filters, e.g., coherence enhancing diffusion
filtering by Weickert (1999). Finally the Hessian
matrix (a 3 3 matrix of partial second derivatives)
has been proposed as a suitable tool in algorithms for
smoothing images of tubular structures, e.g., by Frangi
et al. (1998).
Next, we will describe algorithms that use exactly
these concepts. The difference here lies in the fact that
we will use the results of these algorithms directly
as estimates for the random vector field v(x), rather
than as an intermediate result for subsequent image
smoothing.
Three of the algorithms that will be introduced in
the following sections employ a Gaussian smoothing.
We will choose the known fibre radius r> 0 for the
Gaussian kernel parameter s. This choice is backed
up by results from linear scale space theory Lindeberg
(1994), where it is known that the standard deviation
of the Gaussian kernel is a measure for spatial scale.
Maximumresponseofanisotropic
Gaussians
Let g
u
denote the mask of an anisotropic 3D linear
filter with a prolate shape, u the direction of its longest
prinicipal axis, and  convolution. Then the filter
response (g
u
 f)(x) assumes its maximum when u is
aligned to the preferred fibre direction in a point x2X.
This observation leads to an algorithm for computing
the random vector field v(x): At every point x2X, find
the direction v that maximises the filter response,
v(x)= argmax
u2S
2
+
(g
u
 f)(x): (6)
Since we must maximise over the upper hemisphere
S
2
+
in every pixel, the maximisation in Eq. (6)
can be computationally expensive. Therefore, we
use anisotropic Gaussian filters, for which a fast
implementation of the convolution was obtained by
Lampert and Wirjadi (2006). Anisotropic Gaussian
filters have the form
g
u
(x)=
1
(2p)
3=2
p
detS
u
exp
  1
2
x
T
S
 1
u
x
 ;
for x2R
3
: (7)
Here,S
u
denotes the positive definite 3 3 covariance
matrix determining the shape of the filter mask. To
implement Eq. (6), we use a mask for which all non-
empty level sets have the shape of prolate spheroids.
In this case,S
u
has four degrees of freedom. Different
useful parameterisations ofS
u
exist. We use spherical
polar coordinates to define the fibre direction v(x), thus
the corresponding parameterisation ofS
u
is given.
Let (j;J) denote the spherical polar coordinates
of u, J2 [0;p] the colatitude and j2 [0;2p) the
longitude. Furthermore, let r=s
1
=s
2
> 0 denote the
standard deviation of the Gaussian kernel Eq. (7) in the
two directions orthogonal to u, and lets
3
 s
1
denote
the standard deviation in direction u. Then we get the
following parameterisation for the covariance matrix:
S
u
= R
j;J
0
@
s
2
1
0 0
0 s
2
1
0
0 0 s
2
3
1
A
 R
j;J
 T
;
with
R
j;J
=
0
@
cosj cosJ  sinj cosj sinJ
sinj cosJ cosj sinj sinJ
 sinJ 0 cosJ
1
A
:
For alternative parameterisations of S
u
, see,
e.g., Lampert and Wirjadi (2006); Wirjadi (2009).
In order to implement the maximisation in Eq. (6),
we sample S
2
+
in m points following the numerical
scheme of Fliege and Maier (1999) for S
2
, as it is
readily modified to yield m directions u
i
2 S
2
+
(i =
1;:::;m) free of edge effects. This results in a set of
m= 98 points on the upper hemisphere with nearest
neighbour distances between 12.7
 and 15
 . We refer
to Wirjadi (2009) for details, where also a list of these
vectors can be found.
In summary, an algorithm for estimating the
random vector field v(x) is derived by choosing in
each pixel x the direction v(x) that locally maximises
the result of the convolution of image f with a
set of anisotropic Gaussian filters g
u
i
; i = 1;:::;m.
Of course, the achievable accuracy of this method
depends on the set of pre-computed directions that are
probed, and therefore directly influences the runtime.
Mainaxisofinertia
The basic idea of this method of Altendorf and
Jeulin (2009) is essentially the same as the one
exploited by the just described method: Locally, a fibre
can be approximated by a prolate rotation ellipsoid
whose shorter half-axis length corresponds to the fibre
radius while the direction of the longer axis is the local
fibre direction. Altendorf and Jeulin (2009) derive the
171

WIRJADI O ET AL: The volume-weighted direction distribution
ellipsoid of inertia from directed distance transforms.
The main axis of the ellipsoid then yields the local fibre
direction.
More precisely, from the result of directed distance
transforms in the 26 spatial directions given by the
maximal pixel adjacency, the moments and axes
of inertia are computed. Start with a binarisation
of the fibre system achieved e.g., by a simple
global threshold as described above. Let ‘2 L =
f 13;:::; 1;1;:::;13g index the discrete directions
v
‘
given by the maximal pixel adjacency in a way
that the symmetry v
‘
= v
 ‘
holds. Denote by d
‘
(x)
the distance from x to the next background pixel in
direction v
‘
. That is, d
‘
(x) is the result of the directed
distance transform. Then the sum of the distances
d
‘
(x)+ d
 ‘
(x) yields the length of the chord through
x in direction v
‘
. The moments of inertia M
i jk
(x) are
derived from the endpoints of the centralised chords
P
‘
(x)=
1
2
 d
‘
(x)+ d
 ‘
(x)
 v
‘
by
M
i jk
(x)=
 ‘2L
(P
‘;0
(x))
i
(P
‘;1
(x))
j
(P
‘;2
(x))
k
for i; j;k= 0;1;2;
where P
‘;m
; m= 0;1;2 denote the components of the
P
‘
. Finally, the inertia tensor is given by
M(x)=
1
26
0
@
M
200
(x) M
110
(x) M
101
(x)
M
110
(x) M
020
(x) M
011
(x)
M
101
(x) M
011
(x) M
002
(x)
1
A
:
The eigenvector of M(x) corresponding to its largest
eigenvalue is the main axis of inertia, or an estimate
the local fibre direction v(x) in x.
The direction discretisation emphasizes the 26
used directions and thus results in a systematic error.
Altendorf and Jeulin (2009) derived a correction
reducing both the maximal and the mean error
drastically from about 10
 to less than 5
 and from 6:4
 to 1:3
 , respectively. Introduce the empirical correction
term
t(b)=
(
 0:2sin(p(4b=p)
0:424
); ifb <p=4;
0; otherwise.
:
The corrected direction estimate v
0
is
v
0
= v+
 ‘2L
t(arccos(v
‘
v))kv
‘
 vk: (8)
Note that arccos(v
‘
v) is the angle formed by v and v
‘
andkv
‘
 vk is the length of the orthogonal projection
of v
‘
on the plane perpendicular to v. As mentioned
above, the described version of the algorithm starts
with a binarisation. This can be avoided by utilising
grey-value pseudo-distances as described by Altendorf
and Jeulin (2009).
Thestructuretensor
The fibre direction v is orthogonal to the fibre
surface normals, that is, to the gradient direction in
the image. To account for the circular cross-section
of the fibres, we operate on images that have been
smoothed using an isotropic Gaussian filter. As argued
above, we sets = r, as has also been done by Krause
et al. (2010). That means, the gradient directions at
size scale r are determined. A general problem with
gradients in image data is noise, due to the fact that
they are implemented by difference filters. Therefore
one needs a mechanism for smoothing edge directions.
This can be achieved by considering the tensor product
of gradients, ( f)( f)
T
, a positive semi-definite
matrix with at most one non-zero eigenvalue. The
elements of this matrix are smoothed spatially using
a second isotropic Gaussian filter g
r
:
S
s;r
= g
r
  ( f g
s
) ( f g
s
)
T
 ; (9)
where the outer convolution with g
r
is applied
component-wise, thus resulting in a matrix with
more than one non-zero eigen value, in general. The
resulting S
s;r
is known as the structure tensor and has
successfully been used e.g., in the context of coherence
enhancing diffusion filtering by Weickert (1999). Its
eigendecomposition describes the local image texture.
Especially, the eigenvector with respect to the smallest
eigenvalue of S
s;r
corresponds to the direction with
coherent local grey values, i.e., it points in the direction
where the gradient is minimal. In the context of the
present paper, this is the required local fibre direction
v(x), as the fibre boundaries are expected to produce
large gradients.
From our experience, it is of particular importance
to choose the pre-smoothing parameter s correctly
when applying this algorithm. As it uses gradient
information, it may detect edge structures – instead of
fibres – if the fibre diameter does not fit the chosen
value ofs.
Hessianmatrix
A further method for estimating local fibre
directions is the use of the Hessian matrix. While
gradients contain information on the direction of
edges, second derivatives are known to describe
curvature and therewith the direction of ridge-like
structures (Eberly et al., 1994). Here again, one needs
to handle noise, and therefore we define the Hessian
172

Image Anal Stereol 2016;35:167-179
matrix H
s
computed from a smoothed image f as
H
s
=
  T
 ( f g
s
) (10)
As outlined above, an appropriate choice for the
smoothing parameter is s = r. Using this value,
structures thinner than 2s are suppressed by the
smoothing, see the thinner bar in Fig. 2b. On the other
hand, a much smallers results in erroneous direction
estimates in pixels close to the fibre cores.
Our assumption on the fibre system Eq. (1)
implies in particular that the curvature of the fibres
is small compared to the inverse fibre radius. Thus
the eigenvector of H
s
with respect to the smallest (in
magnitude) eigenvalue is an estimator of the local fibre
direction. Similar observations have previously been
made by Eberly et al. (1994) and Frangi et al. (1998).
(a) (b)
Fig. 2. Demonstration of the effect of smoothing on
the largest (in magnitude) eigenvalue of the Hessian
matrix in a 2D example. (a) Two bars of different
width in some additive Gaussian noise, (b) The largest
(in magnitude) eigenvalue of the Hessian H
s
in every
pixel. The parameter s has been set to correspond to
half the width of the thicker bar.
RESULTS
The applicability of the four different algorithms
for directional analysis of fibre systems is known
from the literature: The Gaussian filter method has
been applied to mCT images of glass and carbon
fibre reinforced composites by Robb et al. (2007),
Wirjadi (2009) and Wirjadi et al. (2009), the main
axis of inertia to glass fibre reinforced composites and
collagen fibres by Altendorf and Jeulin (2009) and
Altendorf et al. (2012), the Hessian matrix was applied
for image smoothing of medical images by Frangi
et al. (1998), and the structure tensor was applied to
images of short fibre-reinforced C/SiC and glass fibre-
reinforced composites by Krause et al. (2010). In the
following, we compare the four algorithms, with the
focus on synthetic data.
RUNTIMES
All four method’s complexity is O(n), where
n is the number of pixels, independent of the
smoothing parameters s and r. For the anisotropic
Gaussian filter method, this is achieved by making
use of the filter separation scheme of Lampert and
Wirjadi (2006). Gaussian filtering can be performed
using O(n) operations using recursive algorithms as
e.g., Young and van Vliet (1995). The 1D distance
transform required as subroutine of the method
of axis of inertia is also known to be of linear
time complexity (Rosenfeld and Pfaltz, 1966). Our
experimental validation of this linear dependency
is based on volumes containing one cylinder with
constant radius (5 pixels) and increasing lengths such
that the volume density is 1%, see Fig. 3.
The results in Fig. 3 confirm the linear complexity,
but show also that the constants ignored inO-notation
differ considerably, here. Especially the anisotropic
Gaussian filter algorithm, which relies on sampling
S
2
+
in m points, is drastically slower than the other
three methods. The Hessian matrix method turns out
to be the fastest one, which can be explained by the
additional spatial smoothing using parameterr that is
applied on all six unique elements of the gradient’s
tensor product to compute the structure tensor, cf.
Eq. (9).
l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l Image size [pixel]
Runtime [min]
l Gaussian
Hessian
Structure Tensor
Main Axis of Inertia
100³ 300³ 350³ 400³ 450³
0 5 10 15
Fig. 3. Runtimes of the four algorithms for image sizes
between 50
3
and 450
3
on an Intel Xeon X5675 at 3.07
GHz. The number of sampling points for the Gaussian
filter method is m= 98, the smoothing parameters are
s = 5 andr = 1 in all cases.
173

WIRJADI O ET AL: The volume-weighted direction distribution
SYSTEMATICERRORSON
SIMULATEDDATA
The usage of simulated image data is the only way
of measuring an image analysis algorithm’s accuracy.
To measure accuracy, we use the difference angle
d = arccos v
T
ˆ v, where ˆ v denotes the mean fibre
direction obtained as the eigenvector corresponding to
the largest eigenvalue of the second order orientation
tensor in Eq. (4).
All image data used in the present section consists
of 128 128 128 pixels, containing one single
cylinder with radius 5 pixels and a length of 80,
centered in the image. The principal direction v of the
cylinder is used as a free parameter in the following
experiments.
We first evaluate the four algorithms’ accuracy
depending on v in order to analyse discretisation
effects. The results are shown in Fig. 4 which
illustrates the systematic errors, where v = (J;j) is
decomposed into its polar coordinate representation
with the colatitude J 2 [0;p] and longitude j2
[0;2p). By keeping one angle fixed and varying the
other, we can asses the methods’ accuracy in this
noise-free setting. Note that strictly speaking, the
results in Fig. 4(a) for J > p=4 lie on the lower
hemisphere. We chose this representation here for an
easier parameterization, nevertheless. By symmetry,
these results (j = p=4 and J > p=4) can easily be
mapped to the upper hemisphere by mirroring the
corresponding vectors on the origin.
The Hessian and the structure tensor as well as the
main axis of inertia methods are capable of measuring
the cylinder’s major axis direction within an error
margin of only a few degrees. When using the method
based on anisotropic Gaussian filtering, on the other
hand, maximal errors are expected when the cylinder’s
direction is far from one of the m samples on the
upper hemisphere. Thus, accuracy of the Gaussian
filter method depends on the number of sampling
directions, i.e., it depends on the discretisation of
S
2
+
in Eq. (6), where a “uniform” distribution of
m points on S
2
+
is a crucial problem. Clearly, with
increasing m, the accuracy of estimation, but also the
runtime increase. This is a clear disadvantage of the
Gaussian filter method. Fig. 4 shows that for some
fibre directions, measurement errors from the Gaussian
filter method even exceed the hemisphere’s maximal
sampling grid distance of 15
 reported by Wirjadi
(2009). It is quite remarkable that the main axis of
inertia method circumvents this trap by the empirical
correction Eq. (8) although being based on just 13
directions.
To give an impression of the algorithms’
robustness with respect to noise, a further experiment
is performed. Straight cylinders model fibres, where
the directions v are chosen independently identically
distributed in S
2
+
. To the binary images of the cylinders,
Gaussian noise with variance s
2
is added. This allows
to plot the resulting accuracyd over the noise variance
s
2
in Fig. 5.
0 50 100 150
0 5 10 15 20 25 30
True theta [deg]
Angular error [deg]
Gaussian
Hessian
Structure Tensor
Main Axis of Inertia
(a) Cylinders parallel to the xz-plane (j = 0). The componentJ of
the axis direction is varied.
0 50 100 150 200 250 300 350
0 5 10 15 20 25 30
True phi [deg]
Angular error [deg]
Gaussian
Hessian
Structure Tensor
Main Axis of Inertia
(b) Cylinders parallel to the xy-plane (J =p=2). The componentj
of the axis direction is varied.
Fig. 4. Systematic error of fibre direction estimation for
cylinders with varying principal directions: Angular
error d of the mean fibre direction computed using
the four considered algorithms from the true, known
fibre direction. The direction of the cylinder’s principal
axis is given in polar coordinates. The smoothing
parameter s is fixed to the cylinder’s radius r, s =
r= 5,r = 1 and m= 98.
174

Image Anal Stereol 2016;35:167-179
0.0 0.5 1.0 1.5 2.0 2.5
Noise variance
Error [deg]
0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90
Gaussian
Hessian
Structure Tensor
Main Axis of Inertia
(a) Angular error of the mean fibre direction as a function of noise
variance s
2
(b) Examples (left: s
2
= 0:5, PSNR=6.0 dB; right: s
2
= 1:5,
PSNR=-3.5 dB)
Fig. 5. Errors of estimated fibre direction from 100
realisations of uniformly distributed directions of
the cylinder under noise. Independent, identically
distributed, additive Gaussian noise with given
variance is used. The fibre radius is fixed to r = 5,
s = r,r = 1 and m= 98.
The results in Fig. 5 show once again that in
noise-free and low noise situations, the Gaussian filter
method’s accuracy is between 10
 and 12
 . For low
noise levels, the main axis of inertia method’s error is
well below 10
 while at a noise variance of 0:75, it
starts to behave worse than the Gaussian filter method.
This is due to the underlying distance measurements
being easily distorted by noise. Switching to the
grey value distance version of the algorithm does not
change the result significantly. The two gradient-based
methods, on the other hand, achieve accuracies of a
few degrees. All methods’ accuracy decreases rapidly
as the noise variance exceeds 1.5 – which corresponds
to a peak signal-to-noise ration (PSNR) of -3.5 dB for
the test data used here. The accuracy of the Hessian
and of the structure tensor methods are the same
(within one standard deviation for most values of s
2
).
In summary, the experiments on synthetic data
presented in this section show that the Hessian matrix
and structure tensor-based estimators are faster than
the main axis of inertia which is in turn faster than the
Gaussian filter-based approach, under all tested fibre
directions v and all tested noise variances s
2
. In terms
of precision, the three methods apart from the Gaussian
filter method performed equally well in our noise-free
experiments, with deviations below 5
 under all tested
fibre directions v. The main axis of inertia is very noise
sensitive. Thus, we finally opt for the faster of the
two grey value derivative based algorithms, i.e., the
Hessian matrix, in the following application examples.
APPLICATION
In this section, we apply the developed methods
to two glass fibre reinforced polymer samples and a
carbon fibre reinforced polymer sample. These two
types of composites which are most common in
industrial applications, today.
Fig. 6. Volume rendering of a (1:5mm)
3
specimen
of glass fibre reinforced polymer, imaged at 3mm
spatial resolution using a laboratory X-ray computed
tomography system.
Fig. 6 shows the volume rendering of a computed
tomography (CT) reconstruction of a glass fibre-
reinforced polymer (GFRP) specimen with 30% glass
fibres(measured by weight). With a known fibre
diameter of 12mm, the second order orientation tensor
of the glass fibre reinforced polymer is estimated from
that data as
175

WIRJADI O ET AL: The volume-weighted direction distribution
A
(2)
=
0
@
0:205  0:020  0:017
 0:020 0:605 0:013
 0:017 0:013 0:190
1
A
: (11)
From the maximum entry in the second diagonal
element of A
(2)
, we can already tell that most
fibres are parallel to y, which we can also compute
via an eigendecomposition of this matrix, yielding
( 0:0505;0:998;0:0335)
T
as mean fibre direction.
While the previously investigated GFRP-sample
was rather homogeneous and the fibres well separated,
the Hessian fibre orientation algorithm will be applied
to a more complex carbon fibre reinforced polymer
(CFRP) sample now. The image data was obtained by
high-quality mCT at beamline ID19 of the ESRF in
Grenoble. The spatial resolution is 0:7mm. In a first
step, we can compute the overall orientation tensor in
the same way as for the GFRP sample above.
Fig. 7. Layers of differently oriented carbon fibres in a
reinforced polymer sample from the aviation industry.
The layers have been segmented in terms of their
angular difference to the global mean fibre direction.
A
(2)
=
0
@
0:200  0:003  0:012
 0:003 0:284 0:020
 0:012 0:020 0:516
1
A
: (12)
Yet, the global tensor in Eq. (12) is not representative
of the entire specimen. If instead of averaging over the
entire image, we pick W in Eq. (4) as a cube with side
length 266mm and compute local tensors in cubes that
cover the image entirely (tiling), we will end up with
different results such as the following two examples:
A
(2)
1
=
0
@
0:207 0:012  0:041
0:012 0:178 0:031
 0:041 0:031 0:614
1
A
;
A
(2)
2
=
0
@
0:198 0:009  0:009
0:009 0:527 0:032
 0:009 0:032 0:276
1
A
:
Of course, the two matrices given above are just
two examples for the results which are available all
throughout the dataset. In order to further visualize
the local differences in the specimen’s fibre system,
we can compare these local orientation tensors to the
global tensor in Eq. (12). We do that by classifying the
entire volume into local areas which are aligned with
the global fibre orientation tensor, and areas which are
not – a process which can be thought of as an image
segmentation algorithm. Using exactly this property,
the layers depicted in Fig. 7 have been segmented. As
a result, we find that the sample is actually composed
of five layers parallel to the yz-plane with differently
oriented carbon fibres: The layers contain fibres which
are rotated roughly 90
 with respect to the fibres in
each neighbouring layer.
As a last application example, Fig. 8a shows the
volume rendering of the CT reconstruction of a GFRP
part from the automobile industry. It was imaged with
a pixel size of 6:7mm and is approximately 6mm wide.
The local fibre orientation analysis yields for the whole
part
A
(2)
=
0
@
0:270  0:016 0:027
 0:016 0:286 0:002
 0:027 0:002 0:444
1
A
: (13)
indicating a strong preference for the z-direction.
However, in the most curved and thin areas, a
significant deviation from this direction can be
observed, see Fig. 8b.
The three application examples that have been
presented here demonstrate the practical applicability
of the methods described in this contribution. They
show that fiber orientation tensors can be useful as a
tool for very detailed, local analyses when inspecting
them individually, as a qualitative tool for visualizing
inhomogeneities such as layers, and as a quantitative
tool when deriving specific quantities such as angular
deviations.
176

Image Anal Stereol 2016;35:167-179
(a) V olume rendering of the fibre component
(b) 2D slice through colour coded deviation from the main fibre
direction
Fig. 8. Part made from glass fibre reinforced polymer.
Volume rendering of the reconstrcuted CT image and
2D slice through 3D map of the deviation from the
main fibre direction( 0:152; 0:026; 0:988). Local
orientation tensors are estimated in cubic subvolumes
of edge length 200mm.
DISCUSSION
We have shown how the volume-weighted
direction distribution computed from an image relates
to the fibre direction distribution of a random fibre
process. This relationship yields estimators not only
for a discretised version of this spherical distribution,
but also for its moments, the so-called orientation
tensors. Furthermore, four numerical schemes for
implementing these estimators have been described.
All of these schemes have the advantage of requiring
only linear time in the number of pixels of a discrete
image. The method based on the Hessian matrix
is quite fast, accurate and robust. It is therefore
recommended to compute volume-weighted direction
distributions.
The limits of the proposed method follow mainly
from the constraining assumption that we have placed
on the dilated fibre processX: Fibres must be solid, that
is neither hollow nor featuring cavities. Fibre cross-
sections have to be circular and of constant diameter.
While the latter could be overcome by applying the
local estimation method for different thicknesses and
finally choosing the “strongest” local fibre direction,
only, the former is essential for all four considered
methods. This assumption will however not be fulfilled
for images of some natural materials such as cellulose
fibres, for instance.
Therefore, the image analysis tools proposed in
the present paper are particularly well suited for
the analysis of glass or carbon fibre reinforced
materials. This follows not only from their geometric
properties. It also follows from the fact that orientation
tensors, which can be directly computed using
the tools described herewith, are important for the
prediction of macroscopic mechanical properties of
such composites.
The comparison clearly favours the two
approaches based on partial derivatives while the
method based on anisotropic Gaussian filters suffers
from the need to sample the unit sphere. It does
however have its merits for 2D applications, see
Schladitz et al. (2016). There, the directional space
corresponds to the unit circle, which is easy to sample
uniformly and runtime limitations are much less of an
issue.
All four methods were applied here “as is”,
that means as described in the respective references.
A hybrid method, combining, e.g., the anisotropic
Gaussian filters with the correction approach used for
the main axis of inertia and the spatial smoothing step
in the derivation of the structure tensor, could further
increase efficiency and robustness.
One of the key assumptions in this presentation is
the constant fibre radius for all fibres. Note that the axis
of least inertia method does not require this. Moreover,
for the partial grey value derivative based methods, it
can be relaxed by varying thes of the Gaussian in the
smoothing step.
177

WIRJADI O ET AL: The volume-weighted direction distribution
ACKNOWLEDGMENTS
This work was supported in part by the
German Federal Ministry of Education and Research
through project 05M2013 (AniS), by the Center for
Mathematical and Computational Modelling (CM)
2
and the Deutsche Forschungsgemeinschaft within the
RTG 1932 “Stochastic Models for Innovations in the
Engineering Sciences”.
REFERENCES
Advani S, Tucker C (1987). The use of tensors to describe
and predict fibre orientation in short fiber composites. J
Rheol 31:751–84.
Altendorf H, Decenci` ere E, Jeulin D, Peixoto P, Deniset-
Besseau A, Angelini E et al (2012). Imaging and 3d
morphological analysis of collagen fibrils. J Microsc
247:161–75.
Altendorf H, Jeulin D (2009). 3d directional mathematical
morphology for analysis of fiber orientations. Image
Anal Stereol 28:143–53.
Altendorf H, Jeulin D (2011). Random-walk-based
stochastic modeling of three-dimensional fiber systems.
Phys Rev E 83:041804.
Bunge H (1993). On-line determination of texture-
dependent materials properties. J Nondestruct Eval
12:3–11.
Clarke A, Archenhold G, Davidson N (1995). A novel
technique for determining the 3D spatial distribution of
glass fibres in polymer composites. Comp Sci Technol
55:75–91.
Eberhardt C, Clarke A (2001). Fibre-orientation
measurements in short-glass-fibre composites. Part
I: automated, high-angular-resolution measurement by
confocal microscopy. Comp Sci Technol 61:1389–400.
Eberly D, Gardner R, Morse B, Pizer S, Scharlach C (1994).
Ridges for image analysis. J Math Imaging Vis 4:353–
73.
Fisher N, Lewis T, Embleton B (1987). Statistical analysis of
spherical data. Cambridge, UK: Cambridge University
Press.
Fliege J, Maier U (1999). The distribution of points on the
sphere and corresponding cubature formulae. IMA J
Numer Anal 19:317–34.
Frangi A, Niessen W, Vincken K, Viergever M (1998).
Multiscale vessel enhancement filtering. In: Proc
Medical ICCA Intervention. Lect Not Comput Sci
1496:130–7.
Hu MK (1962). Visual pattern recognition by moment
invariants. IEEE T Inf Theory 8:179–87.
Kiderlen M, Pfrang A (2005). Algorithms to estimate the
rose of directions of a spatial fibre system. J Microsc
219:50–60.
Krause M, Hausherr J, Burgeth B, Herrmann C, Krenkel
W (2010). Determination of the fibre orientation in
composites using the structure tensor and local X-ray
transform. J Mater Sci 45:888–96.
Lampert C, Wirjadi O (2006). An optimal non-orthogonal
separation of the anisotropic Gaussian convolution
filter. IEEE T Image Process 15:3501–13.
Lee J, Chung M, Alexander A (2006). Evaluation of
anisotropic filters for diffusion tensor imaging. In:
Proc 3rd IEEE Int Symp Biomed Imaging. April 6-9,
Arlington, V A, USA. 77-8.
Lee Y , Lee S, Youn J, Chung K, Kang T (2002).
Characterization of fiber orientation in short fiber
reinforced composites with an image processing
technique. Mater Res Innov 6:65–72.
Lindeberg T (1994). Scale-space theory: a basic tool for
analyzing structures at different scales. J Appl Stat
21:224–70.
Mecke J, Nagel W (1980). Station¨ are r¨ aumliche
Faserprozesse und ihre Schnittzahlrosen. Elektron
Informationsverarb Kyb 16:475–83.
Nagel W (1983). D¨ unne Schnitte von station¨ aren r¨ aumlichen
Faserprozessen. Math Operationsforsch Stat Ser Stat
14:569–76.
Ohser J, Schladitz K (2009). 3D images of materials
structures: processing and analysis. Weinheim: Wiley
VCH.
Ohser J, Steinbach B, Lang C (1998). Efficient texture
analysis of binary images. J Microsc 192:20–8.
Prokop RJ, Reeves AP (1992). A survey of moment-
based techniques for unoccluded object representation
and recognition. CVGIP Graph Models Image Process
54:438–60.
Riplinger M, Spiess M (2012). Asymptotic properties
of the approximate inverse estimator for directional
distributions. Adv Appl Probab 44:954–76.
Robb K, Wirjadi O, Schladitz K (2007). Fiber orientation
estimation from 3D image data: Practical algorithms,
visualization, and interpretation. In: Proc Int Conf
Hybrid Intelligent Systems. Sept 17-19, Kaiserslautern,
Germany. 320–5.
Rosenfeld A, Pfaltz JL (1966). Sequential operations in
digital picture processing. J ACM 13:471–94.
Sandau K, Ohser J (2007). The chord length transform and
the segmentation of crossing fibres. J Microsc 226:43–
53.
Schladitz K, B¨ uter A, Godehardt M, Wirjadi O,
Fleckenstein J, Gerster T et al. (2016). Non-destructive
characterization of fiber orientation in reinforced SMC
178

Image Anal Stereol 2016;35:167-179
as input for simulation based design. Compos Struct To
appear.
Tan J, Elliot J, Clyne T (2006). Analysis of tomography
images of bonded fibre networks to measure
distributions. Adv Eng Mater 8:495–500.
Tankyevych O, Talbot H, Dokladal P, Passat N (2009).
Direction-adaptive grey-level morphology. Application
to 3d vascular brain imaging. In: Proc. 16th IEEE Int
Conf Image Proc (ICIP). November 7-10, Cairo, Egypt
2261–4.
Weickert J (1999). Coherence-enhancing diffusion filtering.
Int Comp Vis 31:111–27.
Wirjadi O (2009). Models and algorithms for image-based
analysis of microstructures. PhD thesis, Technische
Universit¨ at Kaiserslautern.
Wirjadi O, Schladitz K, Rack A, Breuel T (2009).
Applications of anisotropic image filters for computing
2d and 3d-fiber orientations. In: Proc 10th Europ Congr
Stereol Image Anal. June 22-26, Milano, Italy
Yang H, Lindquist B (2000). Three-dimensional image
analysis of fibrous materials. In: Applications of Digital
Image Processing XXIII. July 30, San Diego, USA.
Proc SPIE 4115:275.
Young I, van Vliet L (1995). Recursive implementation of
the Gaussian filter. Signal Process 44:139–51.
Zhu YT, Blumenthal WR, Lowe TC (1997). Determination
of non-symmetric 3-d fiber-orientation distribution and
average fiber length in short-fiber composites. J Compos
Mater 31:1287–301.
179