3D DIRECTIONAL MATHEMATICAL MORPHOLOGY FOR ANALYSIS OF FIBER ORIENTATIONS

In this paper we present algorithms for measuring local char acte istics of random fiber systems. The calculation of the local directions and radii is based on dir ect onal distance transforms and evaluation of the inertia moments and axes of the resulting extremities of the centralized, directed chords. The method provides continuous results while minimizing the runtime by using fe w sampled directions. Furthermore several steps of improvement for the computation of orientation and radiu s information are presented. The algorithms are evaluated using synthetic data and applied to images of real microstructures obtained by computer tomography.


INTRODUCTION
Fiber-reinforced composites are nowadays frequently used for building the enclosure of aircrafts, boats or cars.Our application concentrates on fiberreinforced plastics comprising a polymer matrix reinforced with glass or carbon fibers.The aim of our study is to predict the physical behavior of the material from the knowledge of its microstructure, reconstructed from micro or nano tomography images.One physical property of this very light material is the stiffness, which is highly influenced by the anisotropy of the included fibers.The material will be optimized by changing the parameters of an adapted random geometric model and evaluating the physical properties using numerical simulations.
Local geometric characteristics are essential information for the modeling of random fiber networks.Without this information the virtual material, derived from a random fiber model, is not realistic.It is therefore important to start with a model fitted to the real structure, to adapt the calculation of the physical behavior to the real physical properties of the material and then to start changing the modeling parameters to improve the physical behavior.The most important characteristic for the modeling is the orientation distribution, computed from local information.
The basis of our method is a sufficiently good binarization of a 3D image of a fiber system with solid and not too thin fibers.More precisely, the fiber radius should exceed the length of 3 pixels.Below that thickness, the discretization will have too much influence on the results.On these binarized images (with square or cubic grid), we compute the directional diameters in a fixed amount of orientations (4 in 2D and 13 in 3D) using directional distance transforms.The main inertia axes of the endpoints, given by the local centralized chords, provide an estimate of the local orientation.This estimate is biased towards the sampled directions.We present methods to correct this deviation in 2D and to reduce it in 3D.The inertia moments provide additionally the possibility of estimating the fiber radius and of smoothing the results by making use of the ratio of inertia moments.
Finally, our method is extended to gray level images using thresholded quasi distance.A preprocessing considerably reducing bias in estimation is suggested.
There exist already several methods to compute the local orientation in images, like the Gaussian orientation space by Robb et al. (2007) and the chord length transform by Sandau and Ohser (2007).The chord length transform is not yet studied on 3D images so we will compare our method to the Gaussian orientation space.The main idea of both approaches is to sample a certain amount of directions with different directional operators and referring to the direction with the highest filter response.Thus the results are always limited to the chosen amount of sampled directions and for more exact results their amount has to be increased, resulting in a considerable rise of computation time.Moreover, in order to increase the number of sampled directions, it is necessary to choose a finite number of directions as evenly as possible, which is a nontrivial problem on its own.Our approach avoids these problems as it is fixed to a small amount of sampled orientations (4 in 2D and 13 in 3D), whereas the results are obtained in the continuous Euclidean space.

ANALYSIS ON BINARY IMAGES INERTIA MOMENTS OF DIRECTIONAL DISTANCE TRANSFORMS
In this section the calculation of local characteristics like fiber orientation and radius is treated.The algorithms are based on computing the directional distances to the background for every object point.The sampled directions vs i of the directed distance transform are chosen as the complete neighborhood (in 2D 8 neighbors, in 3D 26 neighbors), see Altendorf (2007).In order to achieve a nearly constant result in a cylindrical fiber, the arithmetic average of the two calculated distances for inverse directions d(vs i ) and d(−vs i ) is considered, which is equal to the half chord lengths defined in Sandau and Ohser (2007) The directed distance transform can be calculated efficiently following an adapted version of the algorithm introduced by Rosenfeld and Pfaltz (1966).The chord lengths are achieved by walking twice through the image: the first time forwards with increasing the distances in the image depending on the distance, assigned to his predecessors in the backward directions (neighbors which have been already visited); the second time by walking backwards through the image, we assign to the pixel the increased distance, assigned to his predecessors in forward directions, if those are inside the object.This algorithm runs in linear time with respect to the number of image pixels and assigns to every foreground pixel the directional thickness of the fiber.
From the endpoints P i = d c (vs i ) • vs i , derived from the centralized distances in all sampled directions, we calculate the moments defined by Duda and Hart (1973).In our case the moments can be reduced to: for 2D and The inertia matrices adapted to our case are: for 2D and for 3D.The inertia moments are the eigenvalues of these matrices and the inertia axes are the eigenvectors as defined by Bakhadyrov and Jafari (1999).Because of the different structure of the inertia matrices in 2D and 3D, the main inertia axis in 2D is the eigenvector to the highest eigenvalue (which indicates the elongation in the according direction), whereas in 3D the main inertia axis is the eigenvector having the lowest eigenvalue (which indicates the inertia by rotating the object around this axis).
The defined main inertia axis gives a first estimate of the fiber orientation, which is however biased towards the sampled directions.

CORRECTING THE BIAS
Evaluation of the presented method shows a certain deviation in the orientation estimate as presented in Fig. 1.By considering the endpoints just in a few sampled directions, those directions receive a high weight.This causes an attraction towards the sampled directions, explaining the deviation.The orientation estimate is perfect in those orientations lying on or in the middle of two sampled directions.This nature of the bias motivated a theoretical study of the problem.The fiber is assumed to be a spherical cylinder with radius r, infinite length and orientation v, represented by the angle θ (in 3D θ and φ , derived from the spherical coordinates).
The centralized distances are given by: from which we can calculate the endpoints To calculate the inertia moments and main inertia axes in 2D, we adapt the existing formula to our case: , (2) Replacing the formula for M pq depending on the fiber parameters r and θ and after applying multiple steps of simplification for trigonometrical functions, it was possible to achieve the following equations, which depend only on the main parameters r and θ : with for i π The maximal possible deviation is limited to 10 • in both cases 2D and 3D.From Eq. 7 it is possible to correct the deviation by inverting it.The orientation can be derived from the estimate θ ′ as follows: The deviation and the corrected angle for the 2D case are illustrated in Fig. 1.
There exist also theoretical solutions for the eigenvalue problem in 3D (Jeulin and Moreaud, 2008).With this, it is possible to deduce an equation for the inertia moments and the inertia vectors.However, we did not achieve a simplification of these complex equations and reduction to the main parameters, which would provide the possibility to correct the orientation.Still there is a way to improve the orientation estimation also in the 3D case.Based on the idea of approximation in 2D, we reduce the deviation of the calculated direction by pushing it away from the closest sampled directions.The attraction to a close sample direction is dependent on the distance between it and the real fiber orientation.Therefore, pushing away the computed orientations from the closest sampled directions is controlled by forces depending on the associated distance.First of all we define the forces, whose formula emerged from several tests based on the two dimensional correction curve.
To complete the approach we need to define the direction in which the force operates.We have defined the force direction to be the projection of the sampled direction on the 2D subspace orthogonal to the calculated orientation v, pn(vs i , v).The approximate orientation v ′ is then calculated as follows: This procedure reduces the maximal error from 9.97 • to 4.78 • and the mean error from 6.40 • to 1.27 • .The reduction of the deviation is visualized in Fig. 2 on the unit sphere in colors from 0 • in blue to 10 • in red.Another aspect which needs to be considered only in the 3D case is the non-consistency of chord lengths for different points in the same fiber.In 2D the chord lengths stay constant for every point inside a fiber, except edge effects at the fiber ends.This situation is shown in Fig. 3(a).In 3D the fiber structure is more complex and therefore the assumption is not fulfilled, as shown in 3(b) for a section of a 3D fiber.We can observe that the chord lengths scale down by moving the point of interest closer to the fiber border.It is thus necessary to find another way to create stable measures for every point in the fiber.Near to the fiber core the directional distances stay stable, that's why we base our calculations not on the measures at the point of interest, but at the center of mass of the extremities, derived from the directional distances at the point of interest.

RADIUS MAPS
The second inertia moment λ 2 ∈ [0.75, 1] can be used in the 2D case to recalculate the fiber radius r based on Eq. 5, There is no equivalent formula in 3D, thus we present a second method to estimate the radius, which can be calculated using the centralized distances d c (vs i ), as they hold already the information of the radius, see Eq. 1: Based on these radius estimates, there are various possibilities to compute the final radius estimate.We have chosen a trimmed mean value: discarding the lower radius estimates reduces wrong estimates due to noise or border regions and discarding the higher radius estimates reduces wrong estimates due to crossing regions.The final estimate is computed as follows: The evaluation showed that this method yields better results than the recalculation from the inertia moments, especially in regions, where fibers cross.

IMPROVEMENT BY SMOOTHING
The resulting direction and radius maps can be smoothed by using a mean filter based on the inertia ratio.The inertia moments indicate the elongation of an object in the direction of the inertia axes.For a ball, all inertia moments are the same, whereas for a fiber the first inertia moment differs significantly from the second (in 3D the second and third are similar).In a point, where two fibers cross, the first and second inertia moments are similar.Therefore we use the ratio of the first two inertia moments to indicate the relevance of the orientation information.The moment ratio for the 2D case (where λ 1 ≥ λ 2 ) is defined as: and for the 3D case (where λ 1 ≤ λ 2 ≤ λ 3 ), we define the moment ratio as: To reduce the difference of orientation and radius information of neighbor pixels we smooth the images using a smoothing filter with a structuring element made of a ball with radius given by the radius map and filter weights given by the moment ratio.It is advisable to apply this smoothing first on the radius map and then on the orientation information to avoid mixing the orientations too much, due to a too large structuring element in crossing regions.

RESULTS
Working on synthetic data and knowing the ground truth, yields the possibility to evaluate the methods with an error histogram.Perfect results would show just one column on 0. The method, which has a higher peak near 0 and decreases faster, provides the better results.The error histogram of the angle maps for 2D and 3D synthetic images (Figs. 4 and 5) shows the improvement between the different steps of our method.Furthermore, we compare our method to the Gaussian orientation space, which uses several elongated Gaussian filters in given directions and assumes the local orientation to be the one, which yields the highest filter response, see Robb et al. (2007).This method can be applied directly to the gray-level images and it is thus not necessary to find a binarization.The results are limited to the chosen directions, whereas our method computes angles in continuous space.That implies that the Gaussian method will need much more directions to achieve comparable results, which increases computation time, especially in 3D.The evaluated error histograms are shown in Fig. 6 for 2D and in Fig. 7 for 3D synthetic images.
On the chosen 3D model in an image of 200 3 pixels, our method finishes in about one minute, whereas the Gaussian method needed two hours for comparable results.

APPLICATION ON DATASETS
In Fig. 8 the method is applied on a SAMimage (Scanning Acoustic Microscopy) of a glassfiber reinforced polymer used for the wheel rim of cars.The sample has a volume fraction of 30% of 1 inch long fibers.Imaged is the projection of a thin slice focussed in a depth of 0.1 mm.
The cutout in Fig. 8d illustrates, that also for very thin fibers we can get a reasonable direction estimate.Nevertheless, as mentioned earlier, in too thin fibers (radius less than 2 pixels) the estimated directions are reduced to the sampled directions.This effect is visible in the direction distribution shown in Fig. 8e.For the 4 sampled directions we get unreasonable high peaks, which are caused by discretization limits.In Fig. 9 we apply our method to the CRP plate.In the direction distribution on the unit sphere, shown in Fig. 9f, the two main distributions from the different layers are indicated by red marks.The θ angle maps can be used to separate the layers.A 3D rendering of the separated layers is shown in Fig. 9e.

ANALYSIS ON GRAY LEVEL IMAGES BY QUASI DISTANCE
In dense parallel fiber networks it appears often, that fibers merge at their borders and that thin frontiers disappear during the binarization process.Applying the directional distance transform directly on the gray level images could yield the advantage to detect these thin frontiers.Three different approaches were chosen.The first possibility to measure distances on gray level images is the quasi distance transform, invented by Beucher (2007), with an additional contrast threshold.The second approach uses Gaussian filters of adapting size.The third approach is based on the comparison of a shape model to the existing fiber structure.

QUASI DISTANCE DEFINITION
The quasi distance evaluates pixelwise for sizes i ∈ N the residual operator τ, derived from the difference between erosions or dilations with a structuring element of varying size i and i + 1 The quasi distance is defined as the size i d for which the dilation or erosion yields the highest residue when compared to the next size i d + 1.In our case, as we want to measure the directional distance, the structuring element is a directed segment.In this case the image can be treated as several 1D signals.For a 1D signal or gray level function f : R + → R we can define the distance for a point x 0 in −X direction with help of the underbuild function This function is increasing and keeps the value in x 0 : fx 0 (x 0 ) = f (x 0 ).The definition of fx 0 equals the reconstruction by dilation from the point x 0 , like it is known in mathematical morphology (see Vincent, 1993 or Salembier andSerra, 1995).The quasi distance is the distance to that point which has the highest gradient: and f c the inverted image (in theory − f , on 8-bit images 255 − f ).
The implementation can be simplified by defining an image walker in the requested direction and buffering the gray level values in a decreasing and in an increasing vector for one line, which needs to be updated respectively.
Furthermore the distance can be influenced by giving a threshold for the significant gradient G 1 .Thus the distance is defined as This threshold treats the case where regions are separated just by a weakly contrasted line and a larger region of background farther is higher contrasted.In the standard case the distance will cross the low contrasted background line and stop at the higher contrasted background.If the threshold is lower than the contrast of the line separating the regions, the distance measure will stop at the line and detect the real fiber end.

GRAY VALUE DISTANCES BY ADAPTING GAUSSIAN FILTERS
This approach makes use of Gaussian filters of adapting size with respect to the distance to the point of interest.Let h(s) = (h i (s)), i = 0, . . ., s be the vector of Gaussian filter weights with σ (s) = (s + 1)/4 and µ = 0.
The filter (with filtersize s x 0 (y) = √ x 0 − y) is applied to the reconstruction by dilation fx 0 with respect to the distance to x 0 : fx 0 (y) = The distance is considered to yield the highest difference in fx 0 (y): with G max g (x 0 ) = sup y∈(0,x 0 ) | fx 0 (y) − fx 0 (y − 1)| .By increasing the size of the filter with increasing distance, the borders are more smoothed in far distances, thus close distances are preferred.

GRAY VALUE DISTANCES BY MODEL COMPARISON
The third approach considers a shape model for the fiber, which takes not only into account the local decrease of the gray level, but also the regularity of the values considered to be fiber foreground.The evaluation of a certain distance h from x 0 is dependent on the regularity of the values between x 0 − h and x 0 and the decrease at the point x 0 − h.The curve is expected to be constantly high on fiber foreground (between x s (h) = (x 0 − h) + s/2 and x 0 ), whereas it should decrease from x s (h) to x e (h) = x s (h) − s.The strength of the decrease can be chosen with respect to the image, on the treated images the minimal choice of s = 2 was optimal.The smoothed model decrease is considered to be like h s (x) = 1 2 sin( (x−(x 0 −h))π s ) + 1 2 .Evaluation is done with x s (h) with the mean value f = x s (h) x 0 f (x)dx, the minimal value the curve does decrease to f min = inf x∈(x e (h),x s (h)) f (x) and the difference between these two values f ∆ = f − f min .The final distance is

EVALUATION ON GRAY VALUE LINES
The presented approaches are evaluated on an original and preprocessed gray value line from the CRP data set.As preprocessing we used toggle mapping to enhance the contrast.The main idea of this filter is to build upper and lower bounds by dilation and erosion, and fit the original gray level in every point to the nearest of the bounds.For more details, see Fabrizio and Marcotegui (2006).Note that this operator can enhance salt and pepper noise.
On the preprocessed values only the thresholded quasi distance and the integral model approach are showing perfect results as presented in Fig. 10.For the approach with the Gaussian filters we can already see that the results are not as stable in neighboring points as for the other approaches.On the original gray value line (without preprocessing) only the approach of the integral model shows acceptable results (Fig. 11).

RESULTS AND COMPARISON
The presented algorithms are tested on a 2D section of a 3D dataset of the carbon fiber reinforced polymer, which was introduced in Fig. 9.During the binarization process, thin contours between fibers are getting lost (see Fig. 12a).Therefore distance measures can cross several fibers, which distorts the calculation results.The thin frontiers between the fibers visible in the gray level image are enhanced by morphological toggle mapping of size 4 (Fig. 12b).For the standard quasi distance it is still possible that the contrast between fiber and thin division line is too low, thus a higher contrast positioned further away is considered as object end.This circumstance causes similar problems in the measurement as in the binary case.The resulting direction maps for the standard quasi distance are shown in Fig. 13b.By using the threshold G ′ = 15 for a sufficient gradient, the thin border lines between the fibers can be detected as object ends and improve the measurements (see Fig. 13c).The comparison of the resulting direction map from Gaussian filters (see Fig. 13d) to the results by thresholded quasi distance is not trivial.The right part seems to have smoother values whereas the left part shows greater deviation to the real orientation.More detailed evaluation can be done from the radius maps presented in Fig. 14.Despite good expectations from the results on the gray value line, the approach of the model shape does not show convincing results in Fig. 13e.The problem with such dense fiber systems is that if the detection of the foreground end fails in just one direction, this direction will carry a too large weight, and intensively influence the estimated orientation.Irregularities in the direction measurement can be smoothed with the adaptive smoothing depending on the moment ratio (introduced in section "Improvement by smoothing".These final results can compete with the result of the Gaussian orientation space applied on the gray level images with filter size (1, 3) (presented for comparison in Fig. 13f. In Fig. 14 the resulting radius maps are presented, which yield the possibility to evaluate the detection of fiber ends.High estimates for the radius indicate errors in the detection of fiber ends.In parallel fiber systems the measurement error caused by merged fibers has a higher influence on the radius calculation than on the direction calculation.Therefore evaluation of the directional foreground end detection is more reasonable on the radius maps.Obviously the approach with the thresholded quasi distance (presented in Fig. 14b) shows the most stable results here.Regarding the computational complexity, the approach of quasi distance which runs in linear time (worst case O(n log n)) yields an advantage compared to the approaches with Gaussian filters (O(n 2 )) or with model shape (O(n 2 log n)).The given times are estimated for one line with length of n pixels.

CONCLUSION
We have seen that the presented method provides stable results, even for the few sampled directions.
The presented analysis tool for binary images works automatically without any parameters and returns maps of local direction and radius estimation.The results are reasonable as shown in the error histograms on synthetic data.Also the computation time is acceptable, for example for an image of 300 3 pixels our algorithm took 9 min, whereas the Gaussian orientation space in fine resolution took about two hours.
The main advantage of the Gaussian orientation space was the direct application on gray level images, for cases where a sufficiently good binarization cannot be achieved.With the thresholded quasi distance method, we have found a reasonable and efficient alternative.Quantitative analysis is in progress.

Fig. 1 .
Fig. 1.Deviation of the direction calculated by inertia moments.True angle of a fiber with length-radius ratio of 1000 vs. calculated and corrected angle.
(a) Deviation of main inertia axes.(b) Deviation of improved orientation.

Fig. 2 .
Fig. 2. Visualization of the deviation of the calculated fiber orientation on the unit sphere.
(a) Constant chord lengths in 2D.(b) Varying chord lengths in section of a 3D fiber.

Fig. 3 .
Fig.3.Illustration of the constant chord lengths in a 2D fiber vs. the variation in a section of a 3D fiber.In both images chord lengths are drawn for two foreground points and the centralized chord lengths are presented additionally outside the fibers, with the main inertia axes in color.

Fig. 4 .
Fig. 4. Error histogram of the angle maps on 2D synthetic data, created by a boolean capsule process.

Fig. 5 .
Fig. 5. Error histogram of the angle maps on 3D synthetic data, created by a cherry-pit cylinder modell.

Fig. 6 .
Fig.6.Comparison of our method (4 sampled directions) to the Gaussian orientation space on a synthetic 2D model, created by a boolean capsule process.

Fig. 7 .
Fig. 7. Comparison of our method (13 sampled directions) to the Gaussian orientation space on a synthetic 3D model, created by a cherry-pit cylinder modell.

Fig. 8 .
Fig. 8. Application to a 2d SAM-image of a glass-fiber reinforced polymer (GRP): (a) initial binary image, (b) the direction map coded using colors, which are explained in (c).(d) zoomed cutout of the image.(e) distribution of the calculated directions.

Fig. 9 .
Fig. 9. Application to a 3D CRP plate (carbon-fiber reinforced polymer): (a) visualization of the whole specimen, where the structure of the four layers is visible, (b) 3D rendering of the treated cutout (cube of 300 pixel side length), (c) structure of the original CT image, (d) colored presentation of the θ angle map (φ is in this material nearly constant): red the background, blue and green the two different layers.(e) rendering of the separated layers, (f) direction distribution on the unit sphere, where red indicates a high and blue a low presence of these directions.
Fig. 10.Comparison of the approaches on a TM4 preprocessed value line.(a) Distances for x 0 = 4 in X-direction, (b) Distances in every point.
Fig. 11.Comparison of the approaches on the original value line.(a) Distances for x 0 = 4 in X-direction, (b) Distances in every point.

Fig. 12 .
Fig. 12. Preprocessed gray level image and binarization of dense parallel fibers in CRP slice.(a) Binary Image, (b) Preprocessed with Toggle Mapping of size 4.