A Monte-Carlo Algorithm for 3D Fibre Detection from Microcomputer Tomography

A model-based approach to analyze fibre distributions in polymer composites applicable for high fibre content is suggested. The algorithm is a four-step iterativemethod usingMonte-Carlo techniques in order to increase speed and robustness for fibre detection. Samples with up to 20% volume fraction of glass fibres and different matrix polymers (PP, PBT) have been analyzed regarding distributions of orientation and length and thickness of the fibres.


Introduction
Short fibre reinforced polymer composites are appealing to mid-sized industry, as their use promises low cost continuous production of parts with enhanced properties. Local fibre length distribution (FLD) and fibre orientation distribution (FOD) [1] determine mechanical properties [2][3][4][5][6][7][8][9] but are strongly affected by production processes [10]. The most commonly used manufacturing process for components made of short fibre reinforced plastics is injection moulding. Hereby, molten plastic material and glass fibres are filled into a mould cavity under high pressure. Particularly in thin walled structural parts, the fibres are oriented in tree layers: two boundary layers and a core layer whose thickness depends on the viscosity of the polymer matrix. The fibres in the boundary layers are mostly oriented in filling direction while the fibres in the core layer are oriented perpendicular to the filling direction. This influences strongly the anisotropy of the composite. Within the development process of moulded plastic parts, flow simulations are performed, for example, in order to avoid weld lines in highly stressed areas by optimization of the positions of the gates. Moreover, the local FOD can also be approximated within such flow simulations. This local FOD information can then be used in anisotropic finite element analysis for the structural part under mechanical loading. However, flow simulations are based on idealized fibre models and lead to uncertainty about FLD as well as FOD [11][12][13]. A comparison of the FOD data computed by flow simulation with real FOD test data is indispensable. Once the test data is available the material parameters, for example, the fibre interaction coefficient in [12], can be calibrated and flow simulation becomes much more predictive. Thus, the usage of short fibre reinforced polymers in safety related components demands effective and accurate quality testing.
While other methods either lack a complete characterization [14,15] or require too much preparation and labour time [16][17][18], microcomputer tomography as a nondestructive testing method is a suitable alternative for those needs [19][20][21][22]. Common fibre contents of commercially available composites vary between 10% and 70% weight. Hence, due to clusters and accumulations, reliable methods will be needed that are able to differentiate and locate single fibres therein.
The system used in this work is a X-ray desktop microcomputer tomography ( CT). Measurements were performed using a SkyScan 1172. Its microfocus sealed X-ray tube operated at 20-100 kV/0-250 A below a maximum power of 10 W. Typical volume picture dimensions used in this work are 1024 × 1024 × 900 cubic voxels of 1.7 m length. This CT together with the algorithm presented in this paper has recently been used in [23] to compute the elastic properties of a short-fibre plastic based on real measured microstructure data by homogenization. In the present work we set our focus on the numerical methodology of fibre detection and discuss the Monte-Carlo algorithm in detail.

Model-Based Algorithm
A common approach of model-based image analysis is to extract image primitives and to establish a correspondence between these and the model primitives. Probabilistic object pattern matching speeds up the scanning process, provides good scaling capabilities if done in parallel, and increases tolerance against noise and distortion. For an overview of image analysis of materials structures, the reader is referred to [24] and the references therein.
The fibre model consists of a chord of cylinders of constant radii and various lengths. The radius of scannable fibres is limited by the minimum cylinder length.
All Monte-Carlo pattern-matching processes run independent of each other. Their results are merged in a controlprocess which detects doublets and optimizes the vector data. Subsequently, this process modifies the voxel data as well as the parameters of the pattern-matching processes and initiates their restart.

Pattern
Matching. Core algorithm steps are Monte-Carlo scanning for separate fibre parts and cylinder integral extrusion for fibre length detection. Scanning for separate fibre parts (index ) is done by calculating radial density functions * ( ) ≡ * in randomly chosen spherical volumes (centres ⃗ ): and is a chosen tolerance.
To approximate radial density, the spherical volume is divided into spherical shell lists of random point-symmetric coordinate pairs within nonequidistant radius-intervals: where denotes the fibre radius and is the pattern radius of the final sphere. Using correlating point-symmetric coordinates according to (2) removes voxels from neighbour fibres as well as noise. Figure 1 shows the radial density * as a function of the integral sphere radius and the cylinder radius , which represents the fibre. The density, that is, the probability that a random point within the sphere is also located within the fibre, is 1 for < . For an increasing sphere radius > we have a decreasing probability and thus a decreasing density distribution which finally tends to zero. The volume of the spherical shell defined by radii and −1 is given by = (4/3) ( 3 − 3 −1 ). And the corresponding volume of the spherical cylinder can be approximated by = 2 ( − −1 ) for > . Thus, the probability that a randomly chosen point can be found within the cylinder shell is * = / and thus which is used as a convergence criterion. Only if this condition is satisfied, a separate fibre has been found and the centre ⃗ of the scanned sphere volume is used for further fibre detection. Otherwise, a spherical volume around another randomly chosen ⃗ within the sphere is checked.
After locating separate fibre parts at ⃗ , the fibre orientation vector ⃗ is calculated either via diagonalizing the matrix of second-order moments or using Jacobian matrices. The fibre length Δ fibre is determined by approximating cylinder integrals using circular coordinate pairs as defined in (6) and using ordinary parameters as angle and radial ( ) and transverse ( ) distances from the origin ( ⃗ ): where the normalized vectors ⃗ , ⃗ 1 , and ⃗ 2 are perpendicular to each other and form a right hand system; that is, For a given starting point and direction, the fibre is made up of small cylinders = −1 + Δ , until the last partial sum Journal of Computational Engineering 3 falls below a certain fraction of the first partial sum Δ < Δ0 ⋅ Δ , which determines the length of the fibre Δ fibre = ∑ Δ . The partial sum is given by After solving the cylinder integrals, all necessary values (position, length, and orientation) for further fibre analysis are available. At this point it is instructive to conclude the complete algorithm for the fibre identification: (1) = 1: chose ⃗ randomly within the CT-image.
(2) Generate concentric shells with radii at the centre ⃗ .
(3) Generate point-symmetric random points within the shells.
(4) Compute mass +1 and centre of gravity ⃗ +1 of all points within the fibre.
(5) Use ⃗ +1 as a new starting point for generating spherical shells.
(7) Compute the radial density distribution of the random points. (4) is fulfilled, a separate fibre is detected. Else = 1 and chose a new ⃗ randomly within the shells and go to (2).
(9) Compute fibre length and orientation by cylinder extrusion.

Duplicate Detection and Match Optimization.
Parallel independent running pattern-matching processes intentionally scan each fibre multiple times starting at random positions. In compounds with high filler fraction, matching quality strongly depends on the starting point of the scan. The control-process has to detect multiple scans of each fibre as duplicates and to select the best match of them. The pattern-matching processes deliver lists of matching information: bottom coordinates ( ⃗ ), directions ( ⃗ ), length (Δ ), radius ( ), and included coloured voxels alias matching hits ( ⃗ 0 , ⃗ ) of (7).
The control-process uses the geometric matching information ( ⃗ , ⃗ , , Δ ) to generate vector representations. Each fibre-match is represented as a pair of nested prisms. The outer prism dimensions are equal to those of the fibre ( ⃗ , ⃗ , , Δ ). Scaling down radius (index ↔) and length (index ↕) of the inner prism is done independently, keeping inner and outer prism-centre upon each other: This prism-pairs are embedded into a 3D-space. To detect duplicate scan-results of a fibre, inner prisms are tested for overlapping. In order to prevent O( 2 ) overlapping tests, all fibres are managed by a binning system, as seen in Figure 2.
The binning system consists of blocks subdividing the test volume. Each binning-cell contains information about all objects touching it. This reduces the number of overlapping tests to the objects touching the same subvolumes.
In case of duplicates, ( ⃗ , ⃗ ) is used as heuristic to sort out the best match.

Iterative Strategy.
The goal of iterative scanning of the volume is to remove well-found fibres, in order to get residual fibres more separated and therefore easier to detect in subsequent scans speeding up the first step of the Monte-Carlo scanning for separate fibres. The procedure is somehow similar to separating bundles in the game of mikado.
As mentioned earlier, the model-based algorithm iterates through a loop of parallel distributed pattern-matching subprocesses, controlled by a centralized main process merging the match-data and modifying the volume image.
The central process also decides which fibres to delete from the volume image and how to adapt the parameters of the distributed pattern-matching before passing them to the subprocesses via network.
In the first iterations a descending fraction, of possible voxels of a matched fibre is used.
In the case | − −1 | < Δ , radius boundaries and length are considered as well.
Depending on how many fibres do contribute to the recent deletion-criteria, parameters for pattern-matching subprocesses are adapted: decrease the pattern radius in order to find separate fibre parts accumulations, increase in (2) in order to tolerate bigger colour-differences, and increase in (6) to enhance tolerance against blurring or cloudy fibres. The number of iterations needed to detect a specific fraction of fibres obviously depends on the image-quality, that is, the contrast between matrix and filler as well on relation  resolution and fibre radius. Figure 3 shows for a voxel size of 1.7 and a fibre radius of 6 that, in case of poly-propylene (PP) containing up to 8%-volume glass fibres, 4 iterations are needed to detect 90% of the fibre volume, whereas polybutylene terephthalate (PBT) compounds with 10%-volume of glass fibres need 7 iterations in order to detect at least 90% of the fibres-volume. Figure 4 shows for PP with 4%-volume glass fibres the time fractions for each step of the pattern-matching of all subprocesses. Calculating the cylinder extrusion integrals x y z Figure 5: Separation of virtually overlapping fibres, Step 1: scanning for fibres (spheres) and direction analysis (blue arrows).
x y z Figure 6: Separation of virtually overlapping fibres, Step 2: calculating cylinder integrals and refining orientation.
is the most time consuming step. The time consumed is proportional to the amount of undetected fibres. The time needed for searching the volume image to find separate fibres is about constant in each iteration.

Algorithm Robustness.
Reconstructed image volume quality depends on density contrast of filler to matrix, filler content, and image resolution. Low density contrast results in blurred object boundaries. High filler content yields to seemingly connected objects. Low image resolution leads to extra or missing voxels on object boundaries. In order to separate seemingly connected objects, each object is scanned multiple times (see Figure 5). This automatically happens due to the fact that all pattern-matching processes are independent and randomly choose their starting points ⃗ .
Accuracy of orientation-and length-detection may vary for each starting point ( Figure 6) and so the controlling process decides upon heuristics regarding which combination of ⃗ 0 , ⃗ , Δ to keep and which one to drop (Figure 7). In order to reduce the effect of blurred object boundaries, cylinder integral directions ⃗ are varied and the maximum of ( ⃗ 0 , ⃗ ) is searched. Using point-symmetric coordinates allows us to separate fibres very efficiently.
In case of low image contrast between filler and matrix, Δ can be reduced to tolerate missing or extra voxels within the cylinder-volume by increasing the acceptance of partial integral sum Δ .

Results
First, the algorithm was tested using computed volume images in order to estimate error magnitudes and effectiveness under various idealized circumstances. After that, short glass fibre reinforced PP composites of 0.5% vol, 1.0% vol, 2.0% vol, 4.0% vol, and 8.0% vol to 14% vol. Finally, short glass fibre reinforced PBT composites including 20% vol were analyzed.

Fibre Detection Quality.
According to the procedure proposed in [25], the algorithm was tested in several steps. It was applied to virtual generated fibre volume images and low filler content compounds and its results were compared to those of manual fibre tracking using slices of tomographic volume images as well as microscopic pictures of residual fibres after removing the matrix material by combustion. Plotting virtual fibres into a volume (see Figure 8) was done to produce a perfectly known system. Additionally this allows us to analyze the effects of different picture defects like blurring and noise against detection efficiency separately.
Compounds with reduced filler content ( Figure 9) were produced to monitor each step of the detection algorithm applied to real fibres. Due to the low filler content, the fibres can be separated by the human eye and this allows manual comparison of detection data to real fibre voxels. In contrast z x y to virtual fibres, real fibre radius may not be constant and their curvatures may not vanish. As shown in Figure 10, aside from pieces of broken fibres, most fibres are detected correctly.
In the next step, 2D-slices ( -layers) of tomographic images of commercial available fibre reinforced composites were used to manually track the fibres in order to compare relative magnitudes of the main components 11 and 22 of the orientation tensor in (11) between the tracking data and the algorithm results.
In Figure 11 results for 11 are shown for several slices along the -axis. In order to check whether fibre length may be an issue for human tracking, the tracked fibres were sorted by length and groups of upper 10%, 20%, 30%, and 50% and all of them were taken into account to calculate the orientation tensor. Within Figure 11, different sizes of the marks "×" have been used to illustrate the different fibre lengths. As can be seen, there is no significant influence on the fibre orientation. In order to assign automatically detected fibre data well defined to -slices, the centre point of each fibre was tested to lie within. To examine scattering of the orientation data, the thickness of the -slices was refined to be 100%, 33%, 20%, and 11% of the gap between the tomographic -layer pictures. Relative thickness z/z max Relative thickness z/z max CT r CT r CT r CT r CT r Figure 11: Comparison of the main orientation tensor element between manual tracking (tomographic slices) and automatic detection of glass fibres (8%-vol) in a poly-amide compound.
Except some slices of the PBT data ( Figure 11), orientation information of tracked fibres and automatically detected fibres are in good agreement.
Manually tracking fibres using microscope images of fibre content after probe combustion provides fibre length statistics to match.

Spatial
Averaging. The fibre vector representations which are generated using geometric matching information ( ⃗ , ⃗ , , Δ ) can be analyzed to get spatial distribution of fibre information ( , , ). One way to achieve this is to divide the test volume into regions. A fibre is taken into account for analysis if its centre lies within a specific region. The region can be any combination of polyhedron and/or spheres. In most cases, the volume is divided into equidistant boxes ( -slices, see Figure 12) in order to examine intermediate layer orientation of fibres within a part.

Orientation Analysis.
In order to describe the fibre orientation state a second-order tensor, = ( is used [26]. Therefore, fibres are characterized as cylinders of equal and constant radius , equal length , starting points ⃗ , and directions ⃗ to determine the components of in (10) as Taking into account variations in radius and length , this yields to a volume weighted form: Figure 13 shows an example of the main orientation tensor components, which can be used directly in this form to validate numerical results from flow simulation.

Length Analysis.
For basic validation of length analysis results, several mixtures of virtual fibres with fixed lengths and correlated random orientation were plotted into a 3Dvolume of voxels. Figure 14 shows an example detection of fibre-lengths with 200, 300, and 400 voxels. In general the peak of the detected lengths lies within an error of 10%, depending on Δ . Overcompensated lowering of Δ due to bad picture quality might yield to erroneous fibre elongation.
In order to validate the resultant lengths on real world probes obtained by the algorithm, we used a crucible to remove matrix material by combustion and viewed the residual fibres under an optical microscope ( Figure 15) where the fibres were tracked manually. Thereby fibres touching the   borders were not taken into account. Figure 16 shows the comparison of length analysis between manual fibre tracking and the computation. As can be seen the length distribution obtained by the proposed algorithm is in a good agreement to the experimental results.

Radius Analysis.
Validating radius results was analogous to the procedure which was used to validate length detections. First virtual fibres were analyzed. Randomly chosen fibres from a small set of configurations (direction, length, and radius) was plotted at different positions inside a voxelvolume. Using the geometric representatives (nested prisms), it is possible to avoid crossings between the fibres. Inverted nested prisms (soft shell inside hard shell) are used in order to adjust minimal radial and tangential distances. Figure 17 shows detected diameter distributions for virtual fibres of 8, 12, and 16 voxels. In general, diameters of perfectly plotted Subsequently, different compounds which were reinforced with the same glass fibre material (original diameter 14 ) were analyzed. Figure 18 shows that the mean diameter generally underestimates the original diameter by half of a voxel size (0.9 ). The mean relative error of the radius detection of fibres approximately corresponds to the picture resolution. For a standard resolution of 1.8 and commercial available standard compounds (diameter 14 ) this results in a relative error of approximate <7%.

Reproducibility
To test the reproducibility of the results for a given tomographic volume image, a poly-amide compound with 8% volume fraction of glass fibres was analyzed 32 times and the probe volume was divided into 15 slices. For each of the 15 slices all elements of the orientation tensor were calculated and an average deviation of all elements for all slices of <0.015 was achieved. The deviation of all deviations was <0.0033.
Analyzing different tomographic volume images of the same specimen again using 15 slices lead to an average standard deviation for all orientation tensor elements of less than 0.025 with a standard deviation of all deviations of less than 0.011.

Conclusion and Outlook
Iterative model-based algorithms to analyze short fibre glass reinforced polymers yield a reproducible and robust method to characterize fibre morphology. The basic algorithm can be adapted to bad picture quality and high filler content, so that partially crossed fibres can be separated directly.
The model-based approach was successfully applied on computed volume graphics as well as on short glass fibre reinforced composites of poly-propylene, poly-amide, and poly-butyl-terephthalate up to commercial filler contents of 30 weight%. In most cases 3 or 4 iterations are sufficient to detect at least 90% of the total fibre volume.
The results of the algorithm were validated by comparing them to manual tracked fibres within tomographic slices and microscopic images. Reproducibility was verified for repetitive analysis of tomographic volume images as well as for different tomographic volume images of the same specimen.
Accuracy and speed of the analysis depend on the aspect ratio of the fibres and on image contrast. Because its statistical basis, the algorithm is highly scalable. Due to its file based protocol, it can be used in heterogeneous clusters already.
In future, the geometric analysis will be extended to support semiautomatic damage examination. Therefore, spatial inhomogeneities of the fibre distribution as well as additional information about transversal and longitudinal distances of the fibres in combination with their length distribution will be used to mark regions of possible fibre fraction. These implementations will additionally be ported to GPUs using CUDA-techniques.