A STOCHASTIC SHAPE AND ORIENTATION MODEL FOR FIBRES WITH AN APPLICATION TO CARBON NANOTUBES

Methods are introduced for analysing the shape and orientation of planar fibres from greyscale images of fibrous systems. The sequence of image processing techniques needed for segmentation of fibres is described. The identified fibres were interpreted as deformed line segments for which two shape and two orientation parameters are estimated by the maximum likelihood method. The methods introduced are shown to perform quite well for simulated systems of deformed line segments with known properties. They were applied to TEM images of carbon nanotubes embedded in polycarbonate.


INTRODUCTION
Fibrous structures are common in biological tissues, like the actin network as a supporting structure in eucaryotic cells, and in industrial materials, like the network of wood fibres in paper.Less known examples include networks of fibres such as, for instance, those formed by electrically conducting fibres, and needed in applications like reinforcement elements, smart clothing, electromagnetic shields or armors, and textile based sensors.Among such fibres the exceptional mechanical and electrical properties of carbon nanotubes make them particularly interesting (Fig. 1).Here they were embedded in a polymer matrix and melt spun to form electrically conducting fibres (Pötschke et al., 2005).
The fibrous shape of carbon nanotubes with very large aspect ratios (length to diameter ratio), even as high as 1000-10,000, allows already at very low volume contents the formation of percolated networks needed for electrical conductivity.Above the percolation threshold the network provides mechanical rigidity and connected pathways for conduction electrons.
The properties of percolating structures of e.g., carbon nanotubes are influenced by their aspect ratio and spatial distribution.In Pegel et al. (2009) clustering of dispersed nanotubes by secondary agglomeration was found to enhance the electrical conductivity of the system.For getting a low percolation threshold, straight and randomly oriented nanotubes are expected to be preferred, but a small degree of orientation in a nanotube system has however been found to be favourable (Du et al., 2005).Even though nanotubes are typically considered as rigid rods, they rather are very flexible, and appear in networks as curved and wavy, which complicates the determination of their orientation.Thus, new methods are needed to properly determine their possible orientation, and thereby the effect of orientation on e.g., the conductivity of the system.
To this end we introduce a stochastic model for the shape and orientation of fibres applied here to carbon nanotubes.In addition, we strive to produce a method for analyzing stochastic fibre systems (Mecke and Stoyan, 1980) in such cases in which the "rose of directions" as a measure of orientational anisotropy does not necessarily provide sufficient means for characterizing the system.
As a stochastic model for individual fibres we used deformed line segments in a similar manner as was used in Grenander and Manbeck (1993) for shape analysis.Further, the collection of fibres was modelled by a Boolean model (Matheron, 1975).To begin with, the model system was composed of a random collection of line segments.For the orientation of line segments we chose a typical model for circular data, i.e., a von Mises density (Mardia and Jupp, 2000) described by two orientation parameters, strength and direction (see an illustration in Fig. 4).When deforming the line segments, each segment was divided into short segments whose orientations were thereafter 'shaked' stochastically such that the distribution of angles between adjacent short segments satisfied a multivariate von Mises distribution with two shape parameters.One parameter described the affinity of a segment angle to that of the previous segment, and the other parameter its affinity to the orientation angle of the original undeformed segment.The resulting 'composite' (deformed) segments were then the fibres of the system and they had stochastically varying shape and orientation (see an illustration in Fig. 5).
A given (two dimensional) image of a fibrous system could then be fitted by the corresponding Boolean model of fibres using the maximum likelihood method by which the parameters of the fibre model, two shape parameters and two orientation parameters, were estimated.This work is an extension to a previous model for individual fibres introduced in Kärkkäinen et al. (2009).
Estimation of these four parameters from the grey-scale data of the image required various image processing phases, e.g., segmentation and identification of individual fibres.In this case the fibres were carbon nanotubes, and we needed to know the pixels that belonged to each nanotube, branching of fibres was not allowed as it is not physically realistic.Therefore, we developed a method for connecting segments of intersecting nanotubes where the resulting structures were physically possible and visually realistic.The method introduced is similar to that in Rizvandi et al. (2008), but we also used the integrated curvature of the fibres and the second derivative of their orientation to properly describe their shape near the intersection areas, in addition to the first derivative already used in Kärkkäinen et al. (2009).
The methods developed were tested against simulated data, and then real data were analyzed, which were taken from carbon nanotubes embedded in a polycarbonate matrix and melt spun to form electrically conducting fibres in microscale (Pötschke et al., 2005).Images of this system were acquired by transmission electron microscopy (TEM) (Fig. 1).For this, ultrathin cuts of about 60 nm thickness were prepared from melt spun fibres of polycarbonate including 2 wt% multiwalled carbon nanotubes, and for the cutting fibres were embedded into an epoxy resin.Cutting was done along the axis of the fibrous material, which was from south-east to north-west in Fig. 1.TEM thus represents a specific 2D projection of the 3D material.Fig. 1.A TEM image of a system of 2 wt% carbon nanotubes embedded in a polycarbonate matrix (the image is reproduced from Pötschke et al., 2005).

IMAGE PROCESSING
TEM images of carbon nanotube systems must be transformed into those of abstract network structures before they could be analyzed by the new segmentation methods developed here.To this end we needed to apply a set of image processing methods, and these are described in this section.
TEM images typically contain imaging noise, and the carbon nanotubes to be identified appeared with an inhomogeneous background, see Fig. 1 for an example.So as to obtain adequate segmentation results, the image was first bandpass filtered (Gonzalez and Woods, 2002).Filtering parameters were chosen so that only details in the size range of carbon nanotube diameter were preserved.To further reduce the noise, a median filter was applied as in Nisslert et al. (2007).This procedure substantially improved the image quality and allowed us to use isodata thresholding in the segmentation of the image into binary form (Ridler and Calvard, 1978).
The binary image was processed by dilation to increase the connectivity of the isolated structures (Gonzalez and Woods, 2002).Since individual small particles do not contribute to the shape of carbon nanotubes, they were removed from the image.After this the image was skeletonized (Gonzalez and Woods, 2002).
The pixels in the skeletonized image were classified into four groups, the background points, skeleton end points, skeleton branch-intersection points and normal skeleton points.This classification was performed using the following ordered rules: If the pixel colour was zero, the pixel was classified as a background pixel.If the 8-neighbourhood of a nonzero-valued pixel contained exactly one pixel with a nonzero value, the pixel was classified as a skeleton end point.If the 8-neighbourhood of a nonzero-valued pixel contained exactly two pixels with nonzero values, the pixel was classified as a normal skeleton point.The remaining nonzero-valued pixels were classified as skeleton intersection points.
After classification each sequence of pixels formed by normal and end points represented a nanotube or part of a nanotube in a bundle of nanotubes.Intersection points were located in the areas where two or more nonparallel nanotubes were touching or one was on top of the other.
We observed that, in the above processing, some physical intersection areas were divided into two or more detached intersection areas.This problem was dealt with such that intersection areas connected with nanotube segments of lengths less than a predefined intersection diameter were merged.The merged areas contained then links to all the remaining segments in the original areas.The short connecting segments were also removed from the image.Furthermore, intersection areas containing no nanotube segments were as well removed.
At this point the actual segmentation of carbon nanotubes could be performed as follows.All intersection areas were processed one by one in a random order.In each case a weight w i j was determined for each possible pair i j of nanotube segments connected to the intersection area.The segments in the pair with the lowest weight were connected to form a single fibre and removed from the intersection.This process was repeated until there were no possible pairs left, i.e., there were one or none segments left in the intersection area.Thereafter the intersection area was considered as segmented and was discarded from further segmentation.
For determination of the above weight, segment directions and shapes near the intersection were quantified.To this end a predefined number of pixels was chosen from each fibre segment starting at the intersection area.Let us mark this set of points for segment i by {p i 1 , p i 2 , . . ., p i n }.For these points, the centroid and principal components were determined.The largest principal component provided a first approximation for the segment-end direction, v i 1 .This set of points was then transformed so that the first principal component became the x-axis, and a parabola y i (x) = a i x 2 + b i x + c i was fitted to the transformed set of points.A second approximation for the segmentend direction was provided by the normalized tangent vector, v i 2 = v i 2 (x), of the parabola at x i 1 such that See Fig. 2 for an illustration of the relevant parameters in segment-end orientation.
Determination of the direction of segment end.Dots mark the centres of adjacent pixels that form an end of a fibre segment.On the top, the first principal component of the set of points is determined and drawn through its centroid.In the middle, the set of points is transformed (T ) so that the first principal component becomes the x-axis.A parabola is fitted to the transformed points.At the bottom, the set of points and the parabola are transformed back (T −1 ).
Note that a parabola could have been fitted also to the nontransformed set of points as a parametric curve.The convergence speed and robustness of the transformation -fit -inverse tranformation sequence was however superior to that without transformation, and the former algorithm was thus much better suited for a completely automatic fitting method.
The shape of segment end was described by the integrated curvature C of the parabola over an interval with The weight for segment pair a, b was then determined such that in which λ 1 , λ 2 and λ 3 are material-specific parameters, C a,b are integrated curvatures of the two intersecting segments and v a,b 1 ( v a,b 2 ) are first (second) approximations to the segment-end directions.The first term on the right hand side of this equation describes the difference in the mean curvatures of the two segments, the second and third terms describe the differences in their directions.The signs of these terms were chosen so that similar curvatures and different directions of the segments resulted in a low weight.This choice favoured the formation of straight fibres.
After identification of the segments connected in the fibre crossings, we could construct sequences of unequally spaced points that could be combined so as to approximate nanotubes in the source image using a linear interpolation technique.For a resulting segmented image, see Fig. 10.
The optimal parameters of the above method (intersection diameter, segment length, λ 1 , λ 2 , λ 3 ) were determined by segmentation of simulated data.To this end the fibre model described in the next section was used to create artificial images for varying fibre coverage, here defined as the percentage of the total (substrate) area covered by fibres.Segmentation result was estimated by comparing segmented simulated data with the corresponding original data with labelled fibres.The number of correctly classified pixels was used as the measure of agreement.Notice that it is possible that the segmentation splits one fibre into two or more regions.In those cases it is ambiguous which of the regions is considered correct.To assure agreement with visual perception, the largest of the regions was chosen as the correct one.Optimization of the agreement with respect to parameters resulted in the optimal parameter vector (intersection diameter, segment length, λ 1 , λ 2 , λ 3 ) = (1, 6, 1.5, 2.5, 1.2).For this set of parameters, the fibre-detection efficiency as a function of coverage is shown in Fig. 3.

MODELLING THE SHAPE AND ORIENTATION OF FIBRES
As a stochastic shape and orientation model for planar fibres we introduced a model of deformed line segments which we call a von Mises fibre.This model is a modification of the shape model used in Grenander and Manbeck (1993).
Let us consider a single line segment that has a fixed length and makes an angle θ 0 ∈ [0, 2π) (the "main angle" in the following) with the horizontal axis.This segment was first 'deformed' such that it was divided into n short segments of equal length, and the orientations of these short segments {θ i } ∈ [0, 2π) were then 'shaked' such that they were distributed by a multivariate von Mises density in the following way.Firstly, the density of angle θ i for a given previous angle θ i−1 and the main angle θ 0 , was given by where the shape parameters, α and β , describe the shape of the fibre such that α measures the 'affinity' of θ i to θ 0 and β its affinity to θ i−1 .In addition, α, β ≥ 0 with α + β > 0. Using standard rules of trigonometric functions and algebra, the right hand side of Eq. 5 can be expressed in the form exp{α cos( Express then (α cos θ 0 + β cos θ i−1 , α sin θ 0 + β sin θ i−1 ) in terms of the polar coordinates where the branches of the arctan function must be taken appropriately.With these definitions the right hand side of Eq. 6 can be expressed as exp[κ i cos(θ i − µ i )] which is an unscaled density of the von Mises distribution.Together with Eqs. 5 and 6 this yields Here the constant I 0 (κ i ) is a zero-order modified Bessel function of the first kind (Mardia and Jupp, 2000).A similar type of density was used in Hughes et al. (2005).
Due to the Markovian way of shaking the orientation angles of subsegments, the density of angles (θ 1 , . . ., θ n ) for a given θ 0 can be expressed in the form This density is called a multivariate von Mises density, and it differs from the one introduced in Grenander and Manbeck (1993) in which modelling was done of closed curves.Also, here θ 0 is regarded as a non-observable random variable in contrast with Kärkkäinen et al. (2009), where it was considered as an unknown parameter.Further, Breckling (1989) defined a von Mises process of angles with longer dependence structure instead of having θ 0 in the model to be estimated.
Consider now a set of von Mises fibres that have random main angles θ j 0 .Without loss of generality we can assume that these θ 0 s follow a von Mises distribution with density where 0 ≤ θ 0 < 2π, 0 ≤ τ < 2π and κ > 0, see Mardia and Jupp (2000) for this and other possible distributions.Parameter τ stands for a "preferred value" of the main angles and κ measures the strength of concentration of the main angles around τ such that small values lead to nearly isotropic fibre systems.
Combining the models, Eqs. 10 and 11, the density of angles of a single von Mises fibre is given by f (θ 1 , . . . ,θ n , θ 0 ; α, β , κ, τ) Regarding θ 0 as an unobserved random variable, we will focus in the following on the density which results from averaging over this angle, in which the integration will be performed numerically.
In order to estimate the parameters of a set of m von Mises fibres described above, we introduce the following notations: Let θ j 0 be the main angle of the jth von Mises fibre and θ j = (θ j 1 , . . ., θ j n j ) the angles of its segments, where n j is the number of these segments.Combining Eqs. 9, 11 and 13, the density of the observable angles (θ 1 , . . . ,θ m ) of all of the von Mises fibres is given by In practice, we can consider any image data of fibrous systems as a realization of a random collection of fibres, which can be modelled using a Boolean model (Matheron, 1975) of deformed line segments as described above.The positions of the original line segments (von Mises fibres) derive from a stationary Poisson point process with an intensity given by the mean number of points in unit area.This means that the number of fibres can be considered as a random variable.The von Mises fibres may also have a random number of segments of equal length, with the density of Eq. 12 for their angles.

ESTIMATION OF PARAMETERS
We estimated the shape (α, β ) and orientation (κ, τ) parameters (generating θ 0 ) from a measured greyscale image by fitting the density of the observable angles by Eq. 14.This estimation required an effective use of the image processing procedures described above for obtaining the needed outlines of single fibres.Furthermore, determination of the division into segments of each outlined fibre of varying length had to be performed.A similar type of approach has been introduced in Grenander and Manbeck (1993) for modelling the shape of potatoes.
After image processing, we had a collection of pixel coordinates for the outlines of m fibres ( j = 1, . . ., m).Using these coordinates, we divided the outline of each fibre into segments as follows.For a chosen end point z j 1 of fibre j, pixel points z j 2 , . . ., z j n j +1 , were searched such that the distance between points z j i and z j i+1 along the fibre was (almost) constant (l) for all i.We thus obtained n j segments.Then the angles θ j i of the segments were determined from the orientations of the lines connecting the end points z j i and z j i+1 of the segments (i = 1, . . ., n j ).For the segment length we chose l = 5 pixels so as to be able to apply the continuous density of angles in [0, 2π) introduced above.
For the estimation of the parameters, we used the maximum likelihood method.The observable angles of the von Mises fibres follow the combined density of Eq. 14.Further, the log-likelihood of the density Eq. 14 is given by l(α, β , κ, τ) containing thus only four parameters that must be estimated.This is a clear advantage over the estimation of likelihood applied in Kärkkäinen et al. (2009).The actual estimation was based on using the R-function optim (R Development Core Team, 2008).

ANALYSIS OF SIMULATED DATA
We can also use the von Mises fibre model to produce simulated data with known properties, which can be used to evaluate the performance of the analysis methods developed.
When producing simulated data, we applied the Boolean model of von Mises fibres as follows.First, a realization of the Boolean model of line segments was produced by numerical simulation.Each line segment had a fixed length of 85 pixels, and the main angles θ 0 of these segments followed the von Mises density Eq. 11 with (τ, κ) = (2.2,10).The locations of the end points (one end) of the line segments, i.e., the starting points of von Mises fibres, were uniformly distributed in the sampling window of size 1024 × 1024 pixels with an edge length of 85 pixels.309 starting points and related line segments were created in that window, and they are shown in Fig. 4. Then von Mises fibres were created from the above line segments.Each line segment was divided into 17 segments of a length of 5 pixels.The angles {θ 1 , . . . ,θ 17 } of these segments were generated from the multivariate von Mises density Eq. 10 utilizing the von Mises densities of Eq. 9.Here we used the shape parameters α = 0.5 and β = 7.0.For an illustration of the resulting set of von Mises fibres see Fig. 5. Since in practice we most often deal with binary or greyscale images, we transformed the contents of the sampling window into a binary image of 1024 × 1024 pixels (not shown).
Assessment of the performance (log-likelihood) of the analysis methods could directly be based on the simulated angles ({θ j i }, i = 1, . . ., 17, j = 1, . . ., 309) of von Mises fibres in the whole area of Fig. 5. Using the log-likelihood of Eq. 15, we thus obtained for the estimated parameters (standard errors were calculated using the inverse of an information matrix in a standard way): α = 0.57 (0.07), β = 7.11 (0.14), κ = 8.60 (0.54) and τ = 2.22 (0.03).It is evident that the parameter values determined through the analysis tools developed are fairly close to their input values.We then considered a binarized version of the simulated set of von Mises fibres.Using the image processing techniques descibed above, fibres of the binary image were segmented from the background, and the outlines of single fibres were determined.The starting point of a fibre was chosen such that was maximized with p the known production direction of the material, i.e., the average direction of fibres; x c is the centroid of fibre pixels and x the chosen starting point of the fibre, which can be at either end of the pixel sequence.As described in the previous section, division into segments was performed starting from the chosen starting point.As the result, we obtained 261 single fibres with a starting point inside the sampling window of   From the determined segments of the fibres (Fig. 7) the angles {θ j i } were determined so as to be able to estimate the parameters.Maximizing the loglikelihood of Eq. 15, we found for their estimates with standard errors: α = 0.42 (0.10), β = 8.49 (0.21), κ = 3.98 (0.42) and τ = 2.24 (0.04).In comparison with the previous (ideal) case, the estimated values are less accurate.This was also expected as now the image was first binarized and then the single fibres together with their starting points were identified.Further, division into segments of the fibres using their outlines is expected to have had a small effect.To this end we made a small simulation study where 2 % of the fibres in Fig. 5 had wrong starting points.We found, as an example, the estimates α = 0.57, β = 7.15, κ = 4.73 and τ = 2.21, from which we can see that κ is sensitive to the selection of starting points.In order to assess the decrease of κ, we further simulated line segments and von Mises fibres using the estimated values α = 0.42, β = 8.49, κ = 3.98 and τ = 2.24 as the parameter values.It is evident that orientation is weaker in Figs. 8 and 9 than in Figs. 4 and 5.

ANALYSIS OF REAL DATA
In this analysis we used the TEM image of a system of carbon nanotubes shown in Fig. 1.The image processing techniques decribed above were applied to segmentation of single fibres in this image, and segmentation results are illustrated in Fig. 10.Then we proceeded to determine the corresponding von Mises model of the identified fibres.The starting point of a single fibre was chosen such that the p of Eq. 16 described the production direction of the material evident in the image: from south-east to north-west.Having determined the starting points, the segments of the fibres and their angles {θ j i } were determined as described above.The resulting identification of 644 von Mises fibres is shown in Fig. 11.Having now the angles {θ j i }, the shape (α, β ) and orientation (κ, τ) parameters could be estimated with standard errors by maximization the log-likelihood of Eq. 15, with the result α = 0, β = 5.10 (0.15), κ = 2.33 (0.17) and τ = 2.28 (0.04).Using these values, we simulated line segments and von Mises fibres shown in Figs. 12 and 13.As the lengths of line segments, we used the mean length of the identified fibres, which was about 30 pixels.Further, each von Mises fibre had 6 short segments of a length of 5 pixels.Orientation seems to be weaker in the simulated data.

DISCUSSION
New and improved methods were introduced for analysis of the shape and orientation of fibres from greyscale images of fibrous systems.Systems of von Mises fibres were generated numerically to estimate the performance of the methods developed, and as a special example a TEM image of a system of carbon nanotubes was analyzed.
Image processing techniques were introduced for handling noisy greyscale images so that they could be properly binarized, and individual fibres could be segmented from them.In segmentation the challenge was to handle fibre intersections, and to this end a rotationally invariant method was introduced.The accuracy obtained in segmentation using these techniques compared well with the present state of the art, see e.g., Rizvandi et al. (2008).Note that the segmentation algorithm introduced does not try to guess fibre paths through intersection areas.It is evident that also these new methods work best for low density systems, and that the easiest way to improve the result of segmentation is to increase the image resolution and to decrease the inhomogeneity of the background.
To analyze the shape and orientation of the segmented fibres, they were interpreted as von Mises fibres.Thereafter their shape and orientation could be estimated based on the maximum likelihood method.The fibre model introduced is an application of the model of deformed line segments.It included two orientation and two shape parameters which were estimated.This estimation was much faster than that in Kärkkäinen et al. (2009), in which a somewhat similar approach was introduced, but the number of parameters was there very large.
The performance of the shape and orientation analysis methods developed could be evaluated using simulated systems on von Mises fibres with known parameters.In the ideal case, when the angles of the simulated segments were directly used, the estimated and input values of the shape and orientation parameters were in good agreement.When various image processing techniques were used before the analysis, the results of estimation were less accurate very much as expected.When analyzing binary or greyscale images, identification of fibres and their starting points, and division into segments of fibres using their outlines, all affect the results.The most sensitive parameter to defects in the analysis appeared to be the one that describes the strength of orientation, while the orientation direction appeared to be rather insensitive to such factors.Both shape parameters behaved in a rather similar manner, and their values could be estimated with a 15-20 percent accuracy even after image processing was applied.
Using simulation studies, we showed that the selection of starting points is critical for the estimation of orientation strength.In practical applications, it will be important to utilize all possible prior information of the system as an input so as to reduce the size of the estimation problem.In the application to carbon nanotube systems considered here, the overall process direction, deduced very easily from the image, made the selection of the starting points of the fibres a much easier and more accurately performed task than otherwise would have been possible.
In the future work, a fibre model less sensitive to the selection of starting points will be sought for.

Fig. 3 .
Fig. 3. Detection efficiency as a function of coverage for simulated data.The coverage in the real carbon nanotube sample is also indicated.

Fig. 4 .
Fig. 4. Simulated line segments with a preferred main angle that indicated orientation of fibres in the north-west direction (orientation parameters (τ, κ) = (2.2,10)).The starting points of the von Mises fibres are denoted by circles.The square is the sampling window of a size of 1024 × 1024 pixels.
Fig. 5.The identified fibres are shown in Fig. 7, and Fig. 6 illustrates how well individual fibres of the image were identified in one of the most challenging locations in the image.

Fig. 6 .
Fig. 6.Part of a binary image of simulated fibres corresponding to the upper-middle part of Fig. 5 (top-left) and the corresponding image of identified fibres coloured with random colors (topright).The same region coloured according to true fibre paths (bottom-left) and the identified fibres with the following markers (bottom-right): circles depict correctly identified fibre intersections, rectangles depict incorrectly identified intersections, black circles depict correctly identified starting points of fibres and black triangles depict incorrectly identified starting points of fibres.

Fig. 7 .
Fig. 7.The single fibres together with their starting points identified from a binary version of an image of simulated von Mises fibres.The fibres are divided into segments of a length of 5 pixels.

Fig. 8 .
Fig. 8. Simulation result using the parameter values (τ, κ) = (2.24,3.98) estimated from the binarized image.The starting points of the von Mises fibres are denoted by circles.The square is the sampling window of 1024 × 1024 pixels.

Fig. 10 .
Fig. 10.Part of Fig. 1 (left) and the corresponding image in which identified fibres appear with random colouring.

Fig. 11 .
Fig. 11.The single fibres together with their starting points identified from the TEM (greyscale) image of carbon nanotubes shown in Fig. 1.

Fig. 12 .
Fig. 12. Simulated line segments using the parameter values (τ, κ) = (2.28,2.33) estimated from the real data.The starting points of the von Mises fibres are denoted by circles.The square is the sampling window of 1024 × 1024 pixels.