Image Anal Stereol 2013;32:13-25 Original Research Paper
AUTOMATIC SEGMENTATION AND STRUCTURAL CHARACTERIZATION OF LOW DENSITY FIBERBOARDS
JÉRÔME LUXB
Laboratoire des Sciences de l'Ingénieur pour l'Environnement (LaSIE) FRE-CNRS 3474, Université de La Rochelle, Avenue Michel Crépeau, 17042 La Rochelle Cedex 1 e-mail: jerome.lux@univ-lr.fr
(Received April 24, 2012; revised February 6, 2013; accepted February 9, 2013) ABSTRACT
In this paper, a new skeleton-based algorithm for the segmentation of individual fiber in 3D tomographic images is described. The proposed method is designed to deal with low-density materials featuring fibers of varied sizes, shapes and tortuosities, like composite fiberboards used for buildings insulation. To this end the paths of the skeleton are first classified according to their connectivity, the connectivity of their adjacent nodes, their orientation, their average radius and the variation of the distance transform along each path. This allows for the identification of spurious paths and paths linking two fibers. Reconstruction of the path of the fibers is done thanks to an optimal pairing algorithm which joins paths that show the most similar orientation and radius at each node/link. The segmented skeleton is finally dilated by means of a growing algorithm ordered by the average radius of the fibers in order to reconstruct each identified fiber. As an application, the algorithm is used to segment a 3D tomographic image of a hemp polymer fiberboard designed for buildings insulation. Information such as the number of contacts, tortuosity, length, average radius, orientation of fibers are measured on both the segmented skeleton and reconstructed image, which allow for a thorough characterization of the fiber network.
Keywords: composite fiberboards, hemp fibers, image segmentation, medial axis.
INTRODUCTION
The macroscopic properties of fibrous materials (e.g., thermal conductivity, permeability, mechanical properties) are related to the nature of the fibers on one hand and to the microstructural characteristic of the fiber network on the other hand. Most important characteristics are the effective density, the number and the nature of fiber-fiber contacts, and the orientation, shape and tortuosity of the fibers. Although a lot of structural parameters are readily accessible without any prior segmentation, e.g., porosity, size distributions, local or global orientations (Lux et al., 2006; Altendorf and Jeulin, 2009), the individualization of each fiber allows for complementary measures (like the number of fiberfiber contact) to be carried out on the segmented image. The segmentation of complex fibrous materials is however still a challenging issue, especially when the fibers cover a wide range of sizes and shapes, as is often the case with vegetable fibers.
Existing segmentation methods are either dedicated to fibers having a solid core (polymer, carbon, glass, etc.), or to hollow fibers, like wood or other vegetable fibers. In the case of hollow fibers, the idea is to segment the inner porosity in order to identify each fiber and/or fibers agglomerate (Bache-Wiig and Henden, 2005; Walther et al., 2006; Tessmann et al.,
2010 and Malmberg et al., 2011). These methods are interesting because they allow to study materials with a wide range of densities. However, these approaches may break down if some fibers are heavily damaged or if the fibrous medium is a composite made up of solid-core and hollow fibers. For such composite materials like those we are studying in this paper, it seems more interesting to fill the inner porosity of the vegetable fibers (Lux et al., 2006; Walther et al., 2006; Malmberg et al., 2011) so that it is possible to use algorithms dedicated to the segmentation of solid-core fibers.
Different approaches exist to deal specifically with solid-core fibers. Eberhardt and Clarke (2002) or Tessmann et al. (2010), propose two methods to compute the center-lines of the fibers directly in gray level images, based on the estimation of local orientations. In Tessmann et al. (2010) the Hessian matrix is computed at each voxel in order to find the center points of the fibers (a subset of the medial axis). Fibers path is then reconstructed thanks to a tracing algorithm that follows the minimum eigenvectors. It was applied successfully on densely packed fiber reinforced polymers featuring fibers with a constant and unique diameter. In Eberhardt and Clarke (2002) the medial axis of the fibers is reconstructed with the help of the local orientation vectors (Creighton et al., 1999). In Yang and Lindquist (2000), Lux et al. (2006) and Tan et al. (2006), the segmentation is performed
on the medial axis of the fibers which is a simpler representation of the actual fibers network. At each node of the skeleton, tests are carried out to form pairs of paths that show the most similar orientation. Individualized fibers can then be reconstructed by geodesic dilation of the segmented skeleton in the initial image (Lux et al., 2006). In this work, we are interested in these kind of methods based on the segmentation of the skeleton, which are well suited to deal with low density fibrous materials.
The skeleton (or medial axis transform) is a compact representation of an object, which preserves its shape (with the concept of end-points or anchor points) and its topology, i.e., the skeleton must have the same number of connected components, tunnels and holes as the original object. A skeleton consists of thin paths which join at nodes. Each voxel of a path has at most two 26-neighbors, whereas the voxels of a node have more than two 26-neighbors. Such simpler representation of an image have been widely used as a support for extracting geometrical features and/or segmentation of various kind of materials (Liang et al., 2000; 0ren and Bakke, 2002; Prodanovic et al., 2006). For a set of fibers, it is assumed that a node is generated when either two fibers cross each other or when a fiber splits in two (as it is sometimes the case for vegetable fibers). Existing skeleton-based segmentation algorithms find and merge the correct pairs of paths at each node in order to reconstruct the real path of the fibers. The pairing algorithms generally simply couple two paths that show the most similar orientation as in Lux et al. (2006), with sometimes a limitation of the angular deviation (Yang and Lindquist, 2000). Note that these methods can only be used for low density materials, i.e., when the contact area between fibers is not too important so that the skeleton transform is assured to produce one path per fiber (which for example would not be the case for dense bundle of fibers). These algorithms have several shortcomings:
First, the fibers have a certain thickness and therefore the contact between two fibers often generates an additional skeleton path. Moreover, when the angle between two touching fibers is small, the digitization can also lead to the apparition of several parasite skeleton paths which are needed to preserve the topology. Fig. 1 shows a typical skeleton configuration arising when two fibers touch. The solution proposed by several authors (Yang and Lindquist, 2000; Lux et al., 2006) is to merge nodes that are too close from each other. This very simple method is of course not well suited to handle materials with fibers of different scales because it relies on a fixed length constraint. Worst, this merging algorithm
may generate a snowball effect, whose consequence is an artificial aggregation of nodes if the length parameter is not well chosen.
Fig. 1. Example of a "ladder-shaped" skeleton generated when two fibers (in black) touch. Blue paths are created in order to preserve the number of 6-connected components of the background.
The second problem is related to the over-detection of end-points. Although a lot of methods for computing the medial axis have been devised in the last two decades, they all show to some extent a high sensitivity to surface noise. The consequence is the occurrence of spurious paths that may be a source of artefacts when the skeleton is used for segmentation purpose. Each non significant path indeed generates an extra node, which may trigger the merging of close nodes as explained above. Moreover, the fiber path may be incorrectly detected at intersections of 3 skeleton paths. The commonly adopted solution is to delete any dead-end path of a length smaller than a user defined threshold (Yang and Lindquist, 2000; Lux et al., 2006; Tan et al., 2006). Again, this is of course not suited to deal with a material featuring a large range of fibers radii, as it is the case of some composite fiberboards.
In this work, we propose a more robust method for both the fibers tracing algorithm and the detection of the special skeleton configurations discussed above in the text, without the need to fix an absolute length constraint. This makes our algorithm more suited to deal with materials made up by fibers of very different sizes. To do so, our algorithm relies partly on the distance transform of the original image, in order both to estimate the radius of any skeleton path or node and to track the variation of the distance transform along the paths.
In the first section, definitions and methods used in this work are detailed. Then, the segmentation algorithm is described and illustrated on a simple image in the second section. Finally, the proposed algorithm is applied to a 3D tomographic image of a hemp polymer composite material.
DEFINITIONS AND METHODS
SKELETON PARTITIONING
In this work, the skeleton S (f) of a binary image f is computed with a 6-subiterations parallel thinning algorithm based on the work of Palagyi and Kuba (1998). It features an additional parallel step to delete the remaining surfaces and a final iterative step in order to obtain a thin skeleton. See Lux (2005) for a complete description of the algorithm. The obtained skeleton is thin in the sense of a line of 1 voxel thickness and, when possible, it is centrally located in the sense of the city-block distance (Borgefors, 1984).
Each voxel of the skeleton is classified according to its number of A-neighbors (Liang et al., 2000). A set of voxels P1 is said to be A -adjacent to another set P2 if it contains at least one voxel which is A-adjacent to a voxel in P2. In this work, we use the following definitions (Yang and Lindquist, 2000):
-	A medial axis path point (PP) is a voxel having less than three A-neighbors. A voxel having exactly one A-neighbor is called an end-point (EP). A voxel with no A-neighbor is called an isolated point (IP),
-	A branch-point (BP) is a voxel having more than two A-neighbors.
-	A path P = {p1,...,pn} is a collection of n A-connected voxels having less than three A -neighbors. We note Z—path a path A -adjacent to Z nodes . Note that the number of adjacent nodes of a path is at most 2.
-	A node N = {n1,..., nn} is a collection of n A -connected branch-points. We note Z-node a node A-adjacent to Z paths.
These definitions allow for a first partitioning of the skeleton into paths and nodes based on local adjacency properties. Additionally, we assume that each path can be classified as one of the three following types:
-	Spurious paths: a spurious path is a non significant 1-path, produced by surface irregularities of the fibers.
-	Center-lines: a fiber center-line is the union of several paths A -connected by common nodes. A single 0-path is the simplest possible center-line.
-	Links: a link is the union of a 2-path joining two fibers center-lines with its two adjacent nodes.
Details about the identification of these types of paths are given in the next sections. Once all the paths are correctly classified, the segmentation consists in identifying each complete center-line.
IDENTIFICATION OF SPURIOUS PATHS
Although a lot of pruning algorithms have been developed for 2D and 3D skeleton (e.g., Svensson and Sanniti di Baja, 2003; Serino et al., 2011), we propose here a method specifically designed to prune the skeleton of fibrous objects. The identification of possible spurious paths is based on their connectivity, length, orientation and medialness (relatively to the center-lines of fibers). Spurious paths must satisfy the four following criteria:
1.	Connectivity: they are 1-paths (i.e., paths with an end-point) connected to a node of order 3.
2.	Length: their length is smaller or approximately equal to the radius of the adjacent fiber.
3.	Orientation: it is likely that these paths will show the greatest angular deviation from the others paths connected to the same node.
4.	Decentering: they are preferentially oriented perpendicularly to the iso-distance lines of the distance transform image d (Rosenfeld and Pfaltz, 1968), or are located partly in flat zones (i.e., zones where the gray level is constant) of distance value equal to 1 (i.e., voxels located at the boundary of the fibers), whereas the center lines of the fibers follow the "crest lines" of the distance transform (if we see the gray-level image as a topographic relief), as illustrated in Fig. 2.
Note that the paths must have a length of more than a few voxels to test for criteria 3 and 4 (typically, two or three pixels at least). If the path is too short, then the criteria 1 and 2 are sufficient to decide whether the path should be deleted or not. In certain configurations, the criterion 3 is not always true. In these cases, the fourth criterion is useful to correctly discriminate between a spurious path and a valid center-line.
PATHS LINKING TWO FIBERS
As already mentioned above in the text, when two fibers touch, it may either generate a node or, if the fibers are sufficiently thick, two nodes connected by a
(a)
(b)
Fig. 2. 2D distance transform (a) and distance transform superimposed with the skeleton (b). The paths linking the two center-lines do not follow the maxima of the distance function.
path. Several paths can also appear due to the effect of digitization near the contact area (Fig. 1). In the existing works, these paths and their two adjacent nodes are merged together if their length is smaller than a fixed length threshold. We follow here the same approach as in the above section, in order to do without a fixed length parameter. The differences with the spurious paths lie in the connectivity and the length constraint:
1.	Connectivity: they are connected to two nodes of order 3 (or sometimes more).
2.	Length: they have a length that is approximately equal to the sum of the local radii of the two adjacent fibers.
3.	Orientation: it is likely that these paths will show the greatest angular deviation from the others paths connected to the same nodes.
4.	Decentering: they are preferentially oriented perpendicularly to the iso-distance curves of the distance transform (Fig. 2).
As it was explained previously, if the path is too short (e.g., its length is lower or equal to 2 voxels), criteria 1 and 2 are sufficient to decide whether the path is a path linking two center-lines or not.
ESTIMATION OF THE LOCAL AND AVERAGE RADIUS
The average radius of any object (path, node) may be estimated thanks to the distance transform d (f) of the initial image f. Due to the fact that the medial axis is centered in the sense of the city-block distance,
we use also city-block distance transform throughout this work. Here, the average radius of any object is computed by taking the arithmetic mean of the values of d for each voxel coordinates in the considered object. The radius may be estimated locally by taking only a subset of voxels in the object. Note that the radius of a node is always a local estimation of the radius of its adjacent paths.
This method gives of course the value of the minor axis length in the case of fibers having a non-circular cross section (like, e.g., vegetable fibers). A more relevant approach would be to use a directed distance transform (Altendorf and Jeulin, 2009) instead of the classical distance transform, at the expense of a greater computation time.
PATH-PAIRING ALGORITHM
In the existing literature, the most probable fiber path at each node/link is computed by searching the pair of paths that minimize their angular deviation. As illustrated in Fig. 3, the orientation is not always sufficient to correctly merge paths adjacent to a node/link. In the present work, we therefore add a new constraint based on the relative radius difference.
The algorithm works as follow for any node Nk:
1. Compute the main orientation v over the first nmax (or the total number of voxels of the path if nmax > lB¡) voxels of each path P adjacent to the node Nk. v¡ is the average of the first nmax unit tangent vectors. Note that it is important to limit the number of voxels used to compute the orientation, because we are only interested in the
local properties of the fiber close to the node or cluster. A fiber may indeed be very tortuous and its orientation may change significantly along its path. nmax should however be large enough to achieve an accurate determination of the orientation of the path. In this work, nmax is taken equal to five times the radius of the current path.
2. For each pair of paths P and Pj adjacent to the same node:
(a)	Compute the angular deviation Faij = 1 — , where aij = arccos (V • vj) is the angle between the two local orientations vi and Vj. Note that the orientation vectors vi and Vj are unit vectors pointing toward the node. Fa therefore tends to zero when the two vectors show an opposite orientation (i.e., the angle aij is close to n radians).
\r-—r-\
(b)	Compute the relative difference FRj = 'r.—rj. between the average radii Ri and R j of each path.
(c) If Faij > F™n and FRij < FR: following weighted sum:
form the
Dj (Nk) = (1 - Wr) Fa¡ + WrFR.
(1)
(d) Finally, the algorithm returns a list containing the values of the computed deviations Dij for each possible pairs of paths i and j adjacent to the considered node.
When Eq. 1 is close to zero, it means that the two paths Pi and Pj have a similar orientation and radius. wr e [0,1] is a weighting factor, which may be used to tune the relative importance of the two terms. This parameter has a noticeable impact on the segmentation results and how to fix an optimal value is discussed later in the text. Note that the choice of the two functions Fa and FR appearing in Eq. 1 for quantifying the angular deviation and radii difference is purely arbitrary. These two functions were chosen for their simplicity, but also because of positive tests results. Maybe more relevant functions could be used to achieve a better quantification of the relative difference between two fibers in term of angular deviation and radii difference. F^11™ and FR"ax are optional parameters prescribing a minimum angular deviation and a maximum relative radius difference. These parameters can be useful to limit the authorized change of the fibers orientation before and after an intersection, as in Lindquist et al. (1999), or to prevent the merging of fibers having very different radii.
Fig. 3. Example of a configuration where two skeleton fibers F1 = {P1 U P2} and F2 = {P3 U P4} with a different radius touch in such a way that the angle between the two paths belonging to the same fiber is n/2. In this case, the correct path of the fiber may not be determined with only directional analysis.
The pairs of paths that show the smaller values for Dij are the most probable fiber paths at the considered node. If the number of paths adjacent to the node is odd, it means either that the fiber splits into two paths (which is the case if one of the paths is a non significant path) or that one fiber touch the end of another fiber.
The weighting parameter wr has a significant impact on the result of the path-pairing algorithm, which is used throughout the segmentation process. It is therefore very important to fix its value carefully. Intuitively, wr should be closer to 1 than to 0, because when the radii of two touching fibers are sufficiently different, the knowledge of the orientation of the paths is irrelevant to decide which paths should be merged together. In order to understand how wr impacts on our path-pairing algorithm and to fix an optimal value, let us consider the very unfavorable configuration depicted in Fig. 3 where two fibers touch and bends with an angle of n/2 . Each fiber center-line F1 and F2 is made up by the union of two skeleton paths: F1 = {P1 UP2} and F2 = {P3 UP4}. In order to reconstruct correctly the path of these two fibers, the deviation D12 and D34 should be smaller than the other deviations. In this case the value of the angle a12 between Pi and P2 is n/2 and is the same as the angle a34 between P3 and P4, so we have also a14 = a23 = n. Furthermore, if we assume that RPl ~ Rp2 = R1 and Rp, ~ Rp4 = R2, the deviations may be expressed as follow (Eq. 1):
Dl2 = D34 = 2 (1 — wr) ,
D14 = D23
: Wr
|R1 - R2I R1 + R2
Since the two other deviations are necessarily larger, they are not interesting for our demonstration.
To ensure that our algorithm merges the right paths together, we seek the value of wr such that D12 < D14, which gives:
wr >
Ri i i Ri + 1
3 ^i-l 3 r1 1
if R1 > 1 Ri
(2)
This simple example shows that, in this very unfavorable configuration (i.e., where Fa = 0 for two paths that belongs to two different fibers) , the path-
pairing algorithm gives correct results when
1—wr
2|Ri-Ril
R1+R2
>
. It means that two paths showing the exact opposite orientations (Fa = 0) won't be merged if the relative radius difference is greater than 1-WWl. As the estimation of the average radii of the paths is not always very accurate (especially for the small fibers, due to the digitization and the use of city-block distance transform), the maximum relative radius difference should be quite significant. As an example, for a maximum relative difference of approximately 50% (i.e., = 1.7), wr = 0.65. The value of the Wr parameter can of course be changed depending on the geometrical characteristics of the studied material.
DETECTION OF NON CENTERED PATHS
As mentioned earlier in the text and illustrated in Fig. 2, the center-lines of fibers follow approximately the crest lines of the distance function, whereas non significant paths or paths linking two fibers are oriented perpendicularly to the iso-distance lines or are partly adjacent to the boundary of the object (where the value of the distance function is equal to 1).
To estimate if the path is medially located in a fiber or not, a simple method consists first in computing the number c (Pi) of iso-distance lines crossed along the path Pi as follow:
If Pi = {x1, ••• , xn} is a path of length lPi = n voxels, where two voxels xj and xj+1 are A-adjacent, then we define the local difference Ad (xj) for xj e
Pi \ xn as:
Ad (xj) = |d (xj-+1) — d (xj) I
(3)
where d (xj ) is the value of the distance transform at
voxel xj.
c is then simply the sum of all d (xj) for xj e Pi \ xn and of the number NfiatZones of zones with constant gray-level of values 1 (this last term is of course only
taken into account for fibers with a radius greater than 2):
C (Pi)= £ Ad (xj)+ NFlatZones . (4)
xj ePi\xn
To compare paths with different lengths, we form the following ratio:
C
c (P ) lp — 1
(5)
If we assume that the radius of a path does not vary too much along its length, a low C value means that the path follows the center-line of the fiber, whereas a high C means that the path crosses a lot of iso-distance lines and is probably either a spurious path, or is a bond between two fibers.
SEGMENTATION ALGORITHM
The aim the of segmentation is to retrieve and label all the center-lines and then to reconstruct the whole fibers from these markers. The proposed algorithm implements the criteria defined previously for the detection and the deletion of the spurious paths and links. Fiber tracing is done with the new path-pairing algorithm. Last, the segmented fibers are reconstructed in the original image by a radius-ordered region growing algorithm using the segmented skeleton as seed points.
It is important to notice that the practical implementation of the criteria needed to identify spurious paths and paths linking two center-lines is not straightforward.
The first difficulty lies in the correct estimation of the local radius when the cross section is not circular. As mentioned earlier in the text, the distance transform only gives us the minor axis value in the case of non-circular cross section. In order to correct this estimation, the values of the radii are multiplied by a corrective factor noted r. To set this corrective factor, one can simply measure the ratio between the major axis and the minor axis of a few fibers cross-sections.
The second question is how to use efficiently the criterion 4 and how to fix a threshold value for the decentering ratio C (Eq. 5). Note that a "high" (i.e., close to 100%) C ratio indicates that a path verifying criteria 1&2 (connectivity and length) is very likely to be either a spurious path or a path linking two other center-lines. However, a "low" (close to 0%) C value does not mean that the path is a center-line. As an example, the distance function along a spurious
path that follows the major axis of an ellipsoid cross-section would be rather constant, and therefore the C ratio will be low. Our idea is here to use criterion 4 only when the path to be tested verifies criteria 1&2 (connectivity and length) but do not verify the criterion 3, i.e., when the path shows an orientation similar to one other path adjacent to the same node. The logic is that the orientation vectors can be sometimes inaccurate (especially when the length of the path is small). In this case, a high C ratio indicates that the path is not a center-line. Of course, contrarily to the other input parameters, the threshold (denoted Cmax) value for C has to be fixed quite arbitrarily. The rule is that it should be "sufficiently" high (intuitively higher than 50%) so that the risk of misclassification remains acceptable.
PRUNING
1. For each 1-path Pi of length lPi connected to a 3-node Nk, do in parallel:
(a) Compute the deviations Dij (Eq. 1) for all pairs of paths adjacent to the node Nk . If the path Pi is in the pair that shows the smaller value for Dij (Nk), then set np (Pi) = 1, otherwise, set np (P)= 0.
(b)	Delete Pi if one of the following condition is met:
-	The length of the path lPi is lower than a threshold length lT.
-	lPi > lT and lPi < T1RNi and np (Pi) = 0, where RNi is the radius of the adjacent node.
-	lp > lT and lp < ri RNi and np (P) = 1 and C (Pi) > Cmax, where RNiis the radius of the adjacent node.
(c)	if Bi is deleted, then merge the node Ni and its two remaining adjacent paths together to create a new center-line with a new label.
The threshold length lT may be adjusted by the user but it has to be greater or equal to 2 voxels, simply because a minimum number of voxels is needed in order estimate the orientation or the C-value of the path.
r1 is the radius corrective factor, which must be greater or equal to i.
Cmax is, as explained above in the text, the parameter fixed by the user to decide from which value of C (Pi) the path Pi (which satisfy criterion 1 & 2) must be deleted.
DELETING CONTACT PATHS
1.	For each 2-path Pi of length lp > lT and lp < n (RNk + RN[) connected to two 3-nodes Nk and Nl, with respective radii RNk and RNl, do in parallel:
(a)	Compute the deviations D (Nk) and D (Nl) (Eq. 1) for all pairs of paths adjacent to the two nodes Nk and Nl. Count the number of time, denoted np (Pi), where the path Pi is in a pair of paths that show the smaller value for D (Nk) and D (Nl).
(b)	Temporarily delete path Pi if np (Pi) = 0 or if np (P) = 1 and C (Pi) > Cmax:
(c)	If Pi is deleted, then merge the two nodes Nj and Nk (which are now only 2-nodes) with their respective adjacent paths and create two fibers center-lines with their own label.
2.	Merge the remaining 2-paths connected to two 3-nodes with theia adjacent nodes if lPi < lT or if lPi < (RNk + RN[) and create a new link.
Note that the paths linking two fibers, identified in step 1, are only temporarily deleted. Indeed, while they are not needed in the fiber tracing step, they are kept in the final segmented skeleton in order to preserve the connectivity of the fibrous network for later measurements.
COMPUTING THE FIBERS PATH
Once all the spurious paths or paths linking two center-lines are deleted, the next step is to find the most probable path for each center-line at the remaining nodes or links of order 3 or more. This is achieved with the path-pairing algorithm detailed earlier in the text. Note that, in order to achieve a better pairing, it is important to process nodes of lowest rank first. When the number of paths is equal to 3, it means either that the fiber splits into two paths or that one fiber touch the end of another fiber. In both cases, the remaining path keeps its own label. The routine ends when all nodes or links are processed.
POST-PROCESSING AND RECONSTRUCTION
To limit the over-segmentation, a last postprocessing routine is applied to suppress the centerlines of small length adjacent to only one other centerline. A path Pi is merged with its adjacent path Pj if lPi < r2RN, where RN is the radius of the node connecting the two fibers and r2 is a user controlled parameter used to tune the intensity of the postprocessing.
(a)	(b)	(c)
Fig. 4. A two dimensional test image of several fibers superimposed with its skeleton (a). Figs. 4b and 4c show the segmented skeleton (nodes and links are in pink) and the final reconstructed image respectively. The parameters used in this example are fixed as follow: r = r2 = 2, Cmax = 0.65, ar = 0.75.
In order to reconstruct the individualized fibers, the segmented skeleton is dilated using a region growing algorithm ordered by the average radius of the centerlines. This algorithm works as follow: at each iteration i, only the voxels x of value I (x) such as R (x) > Rmax — i are dilated by the unit ball (26-connexity), where Rmax is the maximum average radius of all the center-lines and R (x) is the radius of the fiber with label I (x). This works well in the case of fibers with circular cross-sections. It could be improved using directed distance transform (Altendorf and Jeulin, 2009) in the case of fibers with different cross-sections, but it would make the whole process more time-consuming.
To illustrate how our algorithm performs, it is applied on a simple two dimensional image featuring fibers with very different radii, ladder shaped skeleton at contact points, fibers that bends at an angle of nearly 90° and fibers splitting in two (Fig. 4a). Figs. 4b and 4c show the segmented skeleton and the final reconstructed image respectively. Note on the one hand that all the contacts paths (in pink in Fig. 4b) and fibers are correctly identified and on the second hand that the radius-ordered growing algorithm allows for an accurate reconstruction of the fibers.
RESULTS
As an application, we used our algorithm to segment a 3D tomographic image of a composite insulator made up of hemp fibers and thermo-fused polyester fibers that ensure cohesion. The image resolution is 8pm/voxel and the dimensions of the studied volume is 512x512x256 voxels (i.e., approximately 4x4x2 mm3). The studied fiberboard
is highly porous (the porosity is slightly greater than 90%), and features two different kind of fibers:
-	Hemp fibers and shives: both are thicker and generally shorter than polymer fibers and show an important variability in terms of diameter and length. Diameters of hemp fibers are comprised between 15 and 200 pm approximately, whereas hemp shives have diameters of the order of a millimeter. Shives show several inner voids (lumen), which is not the case for hemp fibers.
-	Polymer fibers: these are thin fibers showing a rather constant diameter (about 15-25 pm) and a greater length than wood fibers (greater than 1 mm). They can also be very tortuous, as illustrated in Fig. 6.
X
Fig. 5. Macroscopic view of the hemp/polymer composite fiberboard (polymer fibers in white).
Fig. 6. Visualization of the 3D binary image of a hemp polymer fibrous insulator for building insulation (subvolume of size 384 x 384x192 voxels, 8 pm/voxel).
These fiberboards are manufactured thanks to a non-woven textile process, adapted to such fiber based materials. The fibers are randomly oriented before the mechanical consolidation. This process induces a main fibers orientation along transverse planes (perpendicular to the thickness). Yet, a macroscopic view of the fiberboard (Fig. 5) shows that the fibers are not really all oriented on transverse planes, but can show a deviation angle, up to approximately ±45° with the transverse planes.
The binary image is obtained after a filtering step (opening by reconstruction filter; Salembier and Serra, 1995) followed by a region-growing thresholding. In this simple algorithm, the seed points are voxels whose value is higher or equal to the high threshold value. Voxels are added to the regions if they are connected and if their value is higher than a low threshold value. Before computing the skeleton, the lumen of hemp shives are filled using the simple method described in Lux et al. (2006). As illustrated in the 3D visualization of the image (Fig.6), the fibrous network forms a
very complex and interlaced structure, with fibers of very different sizes. The z-axis is oriented along the thickness of the fiberboard and therefore, xOy planes correspond to transverse planes.
A series of segmentation are carried out for different sets of input parameters detailed in table 1. The reference set is fixed as follow: r1 = 3 (because the cross section of hemp fibers is generally anisotropic) wR = 0.65 and Cmax = 50%. Constant parameters are r2 = 4, Fmin = n/2 and lT = 3 voxels. The impact of a change of each input parameter on a few geometric characteristics measured on the output images is summarized in table 1. These results show that the geometric characteristics and the number of detected fibers do not change significantly for all tested parameters. The number of identified fibers vary from 1118 (i.e., a 7% increase relative to the reference set) when r1 = 2 and 1010 (i.e., a 3% decrease relative to the reference set) when Cmax = 0.3. As it is expected, if r is too small, some links and spurious paths are not detected, hence the increase of the number of identified
fibers. Conversely, if Cmax is too small, some paths are wrongly marked as spurious or as links, hence the decrease of the number of identified fibers. At last, note that the impact of wr is quite moderate, because it plays a role only in certain configurations, where the radius of fibers are sufficiently different. The whole sensitivity study, although not presented here, confirms that our algorithm performs well for a reasonable range of input parameters. In the following, we present the segmentation results obtained with the reference set of input parameters.
Fig. 7. Zoom on a volume with lots of interlaced fibers (128 x 64 x 64). Each fiber is assigned a random color.
Fig. 8. A sub-volume of size 256x256x128 voxels showing the segmented fibers.
Fig. 7 shows a small sub-volume of the final segmented image in a zone containing a lot of interlaced fibers (fibers are assigned random colors). Note that despite the complexity and the high connectivity of the fiber network, all the fibers are correctly identified. A larger volume of 256 x 256 x 128 voxels is shown in Fig. 8 for visualization purpose.
Note however that the hemp shives are generally not correctly reconstructed. Indeed, their rough surface produce a very noisy medial axis which is not easy to prune with the pruning algorithm presented on page 19 because of their very anisotropic cross section. The computed average radius of the spurious paths is in fact very close to the radius of their adjacent centerline (length of the minor axis of the cross section). The presence of these spurious paths of similar radii also perturbs the path-pairing algorithm. A possible solution of this problem would be to modify the method for estimating the radius of the paths, using directed distance transform (Altendorf and Jeulin, 2009).
Length and tortuosity of fibers (computed here as the ratio of length of a center-line to the distance between its two extremities) as well as the number of bounds per fiber are measured on the segmented center-lines. The average radius is measured for each segmented center-line as explained in section 2.4. Fig. 9 shows all the distributions measured on the 1045 identified fibers.
The fiber length distribution (Fig. 9a) shows that about 30% of fibers have a length greater than 1 mm, and should be mostly polymer fibers or shives. Fig. 9b shows the distribution of fibers radius. Note that it is not possible here to identify two distinct distributions for polymer and hemp fibers because lots of hemp fibers have a small radius (60% of the fibers have a radius lower than 25 ^m). However, if we plot the relation between length and radius (Fig. 10a), we can notice that long fibers are also always thin fibers and that thick fibers are always short fibers. All these elements show that the two populations of fibers can be clearly identified in the segmented image using both their length and radius. The average fibers length and radius are 972 ^m and 29 ^m, respectively.
In Fig. 9c it is interesting to observe that while most of the fibers have a tortuosity close to 1, some of them can have a very high tortuosity. The average number of bounds per fiber is 3.7 and while most of the fibers have a number of bounds close to this value, some polymer fibers show a very high connectivity, as it is illustrated in Fig. 9d. If we plot the number of bonds per fiber against the length of fibers (Fig. 10b), we can notice that the number of bounds is quite logically directly related to the length of the fibers.
Finally, the reconstructed image is also used to perform directional analysis. The orientation vectors of the fibers are deduced from the moment of inertia tensors of the segmented fibers (note that the results can be quite different if the orientation is computed on the medial axis only). Fig. 11 shows the number of
Table 1. Impact of input parameters on a few geometric characteristics of the output image.
Input parameters
reference set
ri	3	3	3	2	4	3	3
wr	0.65	0.65	0.65	0.65	0.65	0.5	0.8
Cmax	0.3	0.5	0.7	0.5	0.5	0.5	0.5
Results							
Number of identified fibers	1010	1045	1079	1118	1049	1054	1054
Average length (pm)	964	972	973.7	940	983	965	965
Maximum length (pm)	7612	7722	7722	7722	7722	7928	7928
Average radius (pm)	29.3	29	29.4	29.1	29.5	29.1	29.1
Maximum radius (pm)	274	274	241.3	274	274	274	274
fibers that have an orientation vector (6, $) in spherical notations. Here 6 is the azimuth angle (varying from 0 to 360°) and $ is the polar angle (varying from 0 to 90°) defined such that $ = 90° corresponds to the Oz direction. This graphic clearly emphasizes the anisotropic structure of the fibers network, which is linked to the manufacturing process as explained above in the text. It also confirms the macroscopic evaluation of the fibers orientation: indeed, only a negligible proportion of fibers are oriented along the thickness ($ = 90°) and the majority of the fibers form angles $ varying from -45° to +45° with the xOy plane. This is clearly not an optimal structure from a thermal performance point of view. What is more surprising here is that there are less fibers oriented along the y-axis (6 = 90°and 6 = 270°) than along the x-axis (6 = 0°and 6 = 180°). The transverse planes are therefore not isotropic, as would be expected. It is not clear at the moment if this unexpected anisotropy is due to the material deformation during the sample cutting and preparation or if it can be linked to the manufacturing process.
CONCLUSION
This work is an effort towards improving existing skeleton-based methods for fiber identification in low-density fiberboards. To do so, we make extensive use of the distance transform in order to estimate the radius of each path or node of the skeleton. This allows us to propose new definitions for a more robust identification of spurious paths and links, based on the connectivity of the fibers and nodes, the length, radius and variation of the distance transform along each path, as well as an improved path-pairing algorithm which merges the pairs of paths showing the minimal radius and orientation deviation.
As an application, our segmentation algorithm is applied to an image of hemp polymer fiberboard. The results show that the new path pairing algorithm allows to achieve a correct segmentation of complex and interlaced fibers structure (Figs. 7 and 8), which would not be possible with a classical pairing algorithm based on the minimization of the paths orientation only. Both the segmented skeleton and reconstructed image can thus be used to carry out quantitative analysis of geometrical and topological properties of the fiberboard. In particular, some key factors affecting mechanical, thermal or transport properties, like number of bounds per fiber or orientation distribution are readily accessible.
Several limitations of the current implementation are identified: first, the direct estimation of the radius from the classical distance transform is not well-suited to deal with fibers with very anisotropic cross-sections. A solution would be to use directed distance transform to compute the length of the cross-section major axis, instead of the minor axis. Secondly, this kind of approach based on the skeleton is limited to the study of low-density fiberboard. Indeed, when the contact area between the fibers increases (i.e., when fibers form clusters), the medial axis can be reduced to only one main path instead of several segments corresponding to the center-lines of each fibers. In order to still use the proposed method, a new algorithm for identifying the center-lines must be devised. A possible solution could be to use local orientation vectors for example (Eberhardt and Clarke, 2002; Lux et al., 2006; Altendorf and Jeulin, 2009).
ACKNOWLEDGEMENTS
Many thanks to Prof. Christine Delisée and Jérôme Malvestio for their help on the acquisition of the 3D tomographic data.
Sn 0.25
a 0.20
0.15
Length
(a)
(b)
				
	t ■i			
	1			
				
				
1				
	I			
	L			
0.6 -8
0.4 rt
Tortuosity
(c)
¿F
Number of bounds ] (d)
Fig. 9. Distribution of length (a), radius (b), tortuosity (c), and number of bounds per fiber (d).
(a)
(b)
Fig. 10. Bi-variate plots showing the relation between the length and the radius of fibers (a) and the relation between the length and the number of bounds (b).
Fig. 11. Surface plot whose height indicates the
fraction of fibers with an orientation (d, $ ).
REFERENCES
Altendorf H, Jeulin D (2009). 3D directional mathematical morphology for analysis of fiber orientations. Image Anal Stereol 28:143-53.
Bache-Wiig J, Henden P (2005). Individual fiber segmentation of three-dimensional microtomograms of paper and fiber-reinforced composite materials. MSc Thesis, Norwegian University of Science and Technology.
Borgefors G (1984). Distance transformations in arbitrary dimensions. Comput Vision Graph 27:321-45.
Creighton C, Sutcliffe M, Clyne T (1999). The effect of processing on fibre alignment and compressive strength of carbon fibre-epoxy composites. In: Proc ICCM-12, Paris.
Eberhardt CN, Clarke AR (2002). Automated reconstruction of curvilinear fibres from 3D datasets acquired by x-ray microtomography. J Microsc 206:41-53.
Liang Z, Ioannidis M, Chatzis I (2000). Geometric and topological analysis of three-dimensional porous media: pore space partitioning based on morphological skeletonization. J Colloid Interface Sci 221:13-24.
Lindquist W, Yang H, Oh M, Lee SM, Oh W, Venkatarangan A, Shin H, Prodanovic M (1999). 3DMA, a software package for automated analysis of fiber structure in 3D computed microtomography images.
Lux J (2005). Comportement thermique macroscopique de milieux fibreux réels anisotropes : étude basée
sur l'analyse d'images tridimensionnelles. PhD Thesis, Université Bordeaux 1.
Lux J, Delisée C, Thibault X (2006). 3D characterization of wood based fibrous materials: An application. Image Anal Stereol 25:25-35.
Malmberg F, Lindblad J, Östlund C, Almgren K, Gamstedt E (2011). Measurement of fibre-fibre contact in three-dimensional images of fibrous materials obtained from X-ray synchrotron microtomography. Nucl Instrum Meth A 637:143-8.
0ren PE, Bakke S (2002). Process based reconstruction of sandstones and prediction of transport properties. Transport Porous Med 46:311-43.
Palâgyi K, Kuba A (1998). A 3D 6-subiteration thinning algorithm for extracting medial lines. Pattern Recogn Lett 19:613-27.
Prodanovic M, Lindquist W, Seright R (2006). Porous structure and fluid partitioning in polyethylene cores from 3D x-ray microtomographic imaging. J Colloid Interface Sci 298:282-97.
Rosenfeld A, Pfaltz J (1968). Distance functions in digital pictures. Pattern Recogn 1:33-61.
Salembier P, Serra J (1995). Flat zones filtering, connected operators, and filters by reconstruction. IEEE Trans Image Proc 4:1153-60.
Serino L, Arcelli C, Sanniti Di Baja G (2011). On the computation of the (3,4,5} curve skeleton of 3D objects. Pattern Recogn Lett 32:1406-14.
Svensson S, Sanniti di Baja G (2003). Simplifying curve skeletons in volume images. Computer Vis Image Und 90:242-57.
Tan J, Elliott J, Clyne T (2006). Analysis of tomography images of bonded fibre networks to measure distributions of fibre segment length and fibre orientation. Adv Eng Mater 8:495-500.
Teßmann M, Mohr S, Gayetskyy S, Hassler U, Hanke R, Greiner G (2010). Automatic determination of fiber-length distribution in composite material using 3D CT data. EURASIP J Adv Signal Process 2010:545030.
Walther T, Terzic K, Donath T, Meine H, Beckmann F, Thoemen H (2006). Microstructural analysis of lignocellulosic fiber networks. In: Bonse U, ed., Developments in X-Ray Tomography V. Proc SPIE 6318
Yang H, Lindquist WB (2000). Three-dimensional image analysis of fibrous materials. In: Tescher AG, ed. Applications of Digital Image Processing XXIII. Proc SPIE 4115