Image Based Measurement of Individual Fiber Lengths for Randomly Oriented Short Fiber Composites

Among a wide range of fiber-reinforced composites, those with randomly oriented short fibers, which are also known as random-chopped fiber-reinforced composites (RaFCs), are the most common composites owing to its ease of manufacturing, flexibility of composite shapes, and good material properties, including light weight and high stiffness. These properties of RaFCs are involved with the lengths and distributions of fibers inside the composites. However, inspecting the fiber lengths and distribution remains a challenging problem, particularly when the lengths and locations of individual fibers need to be distinguished using only X-ray transmission images. The main difficulty arises from the variety of fiber widths and their frequent intersections. To address this problem, this paper proposes a comprehensive software system to localize fibers and measure their lengths. Our system is inspired by a previous work for tracing human hair strands. To adopt the previous method for RaFCs, our system extends classic Gabor filter to explore the locally best parameter sets to suit different fiber shapes. With this adaptive filter, we can extract the locations and orientations of local fibers more robustly for RaFCs. Then individual fibers are traced by solving an initial value problem of an ordinary differential equation. To avoid erroneous tracing which typically occurs at intersections, our method traces only the non-intersecting parts of the fibers initially. After that, we connect the fiber segments using the proximity of their endpoints and the orientations. Through experimental validations on different fiber samples, we demonstrate the stability of the fiber tracing and the robustness of the fiber length calculation. Our system works properly even for X-ray radiographic images of heavily tangled fibers in carbon-fiber-reinforced thermoplastic laminates taken by X-ray Talbot–Lau interferometer.


Introduction
Fiber-reinforced composite is a type of composite materials that is reinforced by the fibers dispersed within it, and is one of the most standard materials used in a range of industrial fields including aerospace and automotive industries due to its high stiffness in spite of its light weight. Recently, "random-chopped fiber-reinforced composites" (RaFCs) have received much attention for their ease of manufacturing [5,20], resulting in the low manufacturing costs of composite laminates comprising lightweight components. For comparison, three representative forms of fiber-reinforced composites [9] including RaFCs are shown B Siqi Wang siqi@den.t.u-tokyo.ac.jp 1 School of Engineering, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, Japan in Fig. 1, which illustrates thin layers of composite laminates with different arrangements of dispersed fibers. As can be seen, the fibers dispersed in the composites are in either of continuous or discontinuous forms. According to a textbook [8] on the mechanical properties of composite materials, the composites reinforced by discontinuous short fibers, which are called short fiber-reinforced composites (SRFCs), are inferior in stiffness to those reinforced by continuous fibers, but are more preferable in practice to form composites with complicated surface shapes. In addition, because short fibers can be easily embedded into the liquid matrix resin, SFRCs can be produced by simple injection or compression molding. Compared to SFRCs, which require an additional mechanism to align the fiber directions, the manufacturing process of RaFCs is much simpler. Because of these advantageous properties, RaFCs are currently the most widely used fiber-reinforced composites. The material properties of RaFC, such as stiffness, effective modulus, and the capability of load transfer, are characterized by the distributions of orientations and lengths of fibers. Therefore, there have been a large number of studies conducted [2,7,10,21,29] to examine these distributions. For such studies, the measurement of individual fiber lengths is a fundamental task to calculate the distribution. However, in cases where a huge number of fibers with different lengths are involved within a small part, the measurement of individual fiber lengths is still challenging, and a solution for this problem has not yet been well established. One exception was proposed by Terada et al. [29], who leveraged an image obtained by a document scanner to measure individual fiber lengths; however, its applicability to X-ray imaged fibers is unclear because it is tested only with fibers visible on the surface. Lionetto et al. [15] burned the matrix part out from fiber-reinforced composites to produce a sample with only fibers [15], and measured the lengths of individual fibers within images obtained via optical microscopy. de Pascalis et al. [4] proposed a method to extract fibers from three-dimensional (3D) X-ray computed tomography (CT) images using image segmentation methods. However, acquiring a 3D CT volumetric image takes much more time and cost than two dimensional (2D) X-ray radiography; therefore, many studies have only used 2D X-ray radiography to inspect the fibers [25]. For example, such inspections with 2D X-ray radiography have been applied to many RaFCs such as glass fiber reinforced plastics (GFRP) [1] and steel fiber reinforced concretes [11]. Although carbon fibers cannot be imaged clearly by simple X-ray absorption imaging, visibility imaging with X-ray Talbot-Lau interferometer (XTLI) has enabled inspection of carbon-fiber-reinforced thermoplastics (CFRTPs) [18].
Despite the long-standing effort described above, quantitative analysis to calculate the fiber length distribution (FLD) from 2D images has remained a challenging problem due to the difficulty of measuring the lengths of individual random fibers. We review several previous studies that are closely related to ours, although there are a huge number of previous studies on image-based material inspection. The lengths of individual fibers are commonly measured by tracing fiber skeletons, which are given by thinning fiber regions [24] and medial axis transform [32]. The fiber skeletons are traced straightforwardly using the adjacency of pixels in conventional approaches [24,32], graph-based path finding such as the A * method and active contour model, are also applied [31]. However, even using such sophisticated approaches, extracting and tracing dense, intersecting, and randomly oriented fibers is a difficult task. Particularly, the ambiguous orientations and curvatures of fibers at the intersection often cause tracing failure. The connectivity of the fibers at the intersection has been solved in two typical ways in previous approaches. One approach is to connect the fiber during the fiber tracing [32]. When the tracing reaches a junction or an intersection, the currently traced fiber is connected to the one with the most similar orientation or curvature at the endpoint. The other typical approach to resolve the connectivity is fiber clustering [19,24]. In this approach, potential fibers are traced entirely to yield the fiber segments, and then the fiber segments are clustered to form the final fibers using their orientations and curvatures. Unfortunately, both of these approaches are tested in composites with relatively sparse fibers with a few intersections, and their applicability to RaFCs with dense and frequently intersecting fibers must be restricted.
To address this problem, this paper proposes a comprehensive software system for acquiring the fiber lengths of RaFCs by individual fiber measurement. The proposed system is largely inspired by previous studies [3,22] on tracing human hair strands, and is based on the standard techniques of fiber extraction and measurement. However, we make improvements to (i) the estimation of fiber locations and orientations and (ii) fiber tracing of the previous method. In addition, we introduce an additional process of (iii) agglomerative clustering of fiber segments obtained by the preceding fiber tracing process. Figure 2 shows an input image and intermediate results for the three techniques improved by our method. The input to our system is an X-ray radiographic image of a fiber-reinforced composite (see Fig. 2a). The locations and orientations of the fibers are extracted using the standard Gabor filter. In this part, we optimize the parameters of the Gabor filter at each pixel location individually rather than optimizing it globally, which helps to extract fibers with different widths. The locations of fibers are highlighted by the colors corresponding to their orientations in Fig. 2b. Then, individual fibers are traced using the map of fiber orienta- Fig. 2 Input radiographic image and intermediate results of our system for fiber length measurement. Following the standard approach to fiber tracing, fiber orientations are calculated at the beginning, and then the confident parts of fibers are traced to obtain the fiber segments. Finally, the fiber segments are clustered together to form the individual fibers tions. The fiber tracing in our system is formulated as an initial value problem of an ordinary differential equation (ODE), and is solved by numerical solvers such as Runge-Kutta method. The fiber segments obtained by this tracing operation are shown in Fig. 2c. Finally, we perform agglomerative clustering to the traced fiber segments to obtain the final fibers for calculating the lengths. This fiber clustering approach has been used in a previous study [24], wherein we carefully designed a set of criteria to connect fibers at the intersections to achieve sufficiently accurate clustering for dense and random fibers. The fiber segments are clustered to form the final output fibers, as shown in Fig. 2d. To demonstrate the effectiveness of the proposed method, we evaluated our method quantitatively using three handmade examples with the ground truth lengths. We also qualitatively evaluated the results using four real images of CFRTPs taken by using XTLI. As we will show in the experimental results, our system successfully traced and obtained plausible fiber lengths for the dense and intersecting fibers in RaFCs.

Confidence-Based Fiber Extraction
In our system, the fiber locations and orientations are detected from an X-ray radiographic image with "confidences" to describe the reliability of the extracted locations and orientations, which are useful in the subsequent processes of calculating the fiber lengths. The concept of confidence was originally proposed by Paris et al. [22] and more formally defined by Chai et al. [3] for detecting orientations of human hair strands based on the variance of the Gabor filter responses; therefore, the concept itself is not new. However, the fibers of RaFCs have different properties from those of hair strands. We extended this previous work to more robustly detect the properties of fibers, including their orientations and widths. In this section, we begin with the classic problem of detecting fiber structures using 2D digital image processing, and then explain the application and extension of the confidence for more robust fiber extraction.
In Sects. 2 and 3, we use a simpler handmade sample, which is made with fishing lines and Styrofoam, as shown in Fig. 8 to more clearly demonstrate the effects of each computational step in our system. This handmade sample is detailed later in Sect. 4.

Background: Confidence of Fiber Extraction
Detection of thin fiber structures within 2D digital image is a classic problem of the image processing. For example, steerable pyramid [6,28] and Gabor filter bank [12] are the two representative approaches for detecting fiber structures. These fiber extraction methods are traditionally used to obtain local visual features for the image segments applied for texture recognition and segmentation. Therefore, the accurate extraction of individual fibers from digital images is out of the scope of these approaches. Apart from the non-destructive testing, a recent trend in a field of human digitization has considered the importance of extracting individual fibers for the purpose of manipulating and digitizing human hairs [3,16,22]. These approaches used the variance of Gabor filter responses to define how confident the extracted fiber orientations are. To explain these approaches, let us start with the standard Gabor filter kernel.
where u and v are the local coordinates at the filtered domain, andũ = u cos θ − v sin θ andṽ = u sin θ + v cos θ are the local coordinates rotated by angle θ . Further, σ and λ are the parameters to determine the blurriness of the Gaussian filter and the frequency of interest of the cosine part, respectively. Applying this filter at each pixel (x, y) with different θ yields the most reliable orientation of the local structure, which we denote asθ(x, y).
where I denotes the target image to be filtered, * represents the convolution operator, and F(x, y, θ) = |K θ * I |(x, y) is the response of the Gabor filter kernal at pixel (x, y). Then, the variance of Gabor filter responses [22] are defined as follows: The standard deviation of the filter responses are defined ordinarily by a square root of the variance.
This standard deviation is used as the confidence of the detected orientation in the literature [3].

Preprocess
Our fiber extraction starts with filtering an input X-ray radiographic image to enhance the fiber pixels. Typically, our primary target, i.e., the X-ray radiographic image, is more likely to suffer from noises than the portraits captured by ordinary digital single lens reflex (DSLR) cameras used for the hair strands [22]. In contrast, smoothing filters such as Gaussian filter, bilateral filter [30], and non-local means filter [26] can overly blur important details of thin fibers in X-ray radiographic images. Therefore, we instead apply adaptive histogram equalization [33] to adaptively enhance the contrast of fiber pixels. This process is an extension of traditional histogram equalization to enhance the image contrast but preforms the equalization at every local set of pixels. The input image before and after the contrast enhancement are compared in Fig. 3.

Extension of Fiber Confidence
In previous studies [3,22], the confidence of the detected fibers was defined by the variance of Gabor filter responses with respect to only the fiber orientations. However, to inspect the fiber widths and arrangements as well as the fiber orientations, the variance should be calculated with respect to all of these parameters. Based on this concept, we extend the definition of the variance in Eq. (3) by considering also the parameters σ and λ used to define the Gabor filter kernel in (1). A set of best parameters {θ,σ ,λ} for the local region around each pixel is defined to provide the strongest Gabor filter response.
Then, using the best set of parametersΘ, we also extend the definition of the variance in (3) as follows: Here, d angle (θ, θ) is given by Eq. (4). As in [3], we use the standard deviation of the responses, i.e., the square root of the variance in Eq. (7), to define the confidence w(Θ(x, y)) of the parameter set. For simplicity, we hereinafter denote w(Θ(x, y)) as w(x, y). The maximization problem in (6) is a complicated nonlinear and non-convex problem and is difficult to solve using ordinary non-linear solvers, such as Newton method. In this research, we therefore use a simple brute-force approach to search the optimal set of parameters over the discretized parameter combinations. Although the original possible domains of the parameters are θ ∈ [0, π), σ ∈ (0, ∞], and λ ∈ (0, ∞], we narrow the domains of the latter two to be finite for ease of computation. Then, the discrete values are sampled for each parameter. In our implementation, 72 angles are uniformly sampled for θ , and 11 values including both ends of the domains are uniformly sampled for σ and λ. Finally, the best parameter set to maximize (6) is chosen from 72 × 11 × 11 = 8, 712 candidate parameter sets. The domains and discretized candidate values of these parameters for different experimental samples are shown in Sect. 4. In contrast that the fiber widths are thin and almost uniform in the result of the previous method [3], our method succeeds in marking fiber regions with different widths. Further, the regions with confident fiber orientations can be thickened by the diffusion process Eventually, the fiber regions and the corresponding orientations are extracted robustly for fibers with different widths. As shown in Fig. 4, our method successfully marked the fiber regions with different widths in contrast that all fibers are marked with almost the same widths by the previous method by Chai et al. [3].

Diffusion for Fiber Orientation Map
Even though the fibers with different widths are extracted adequately, the orientations of fibers are reliable only at the center of fibers. In other words, regions of reliable orientations are too thin to cause failures in the later fiber tracing process. To alleviate this problem, we "thicken" the fiber regions of the orientation map by diffusing the orientations to the neighbor pixels. For the extracted fiber orientationsθ which ranges from 0 to π , we first compute a structural tensor M for each pixel (x, y).
We then diffuse the four elements of the structural tensor M individually by applying a Gaussian filter to tensor elements weighted by the confidence w.
where ρ denotes a parameter to control the blurriness. Then, we calculate the diffused orientationθ at each pixel by the eigenvalue decomposition of M . According to the definition of M in (8)

Fiber Tracing, Clustering, and Length Measurement
To measure the lengths of individual fibers precisely, our method explicitly acquires the curved shapes of individual fibers by tracing them one by one. The tracing operation of each fiber is equivalent to solving an initial value problem, i.e., an ODE with a set of initial conditions. The ODE is defined by the orientation field, namely the derivative of the function for a parametric curve of each fiber. Let P(t) = (x(t), y(t)) be a function of a parametric curve, where x(t) and y(t) are the coordinates of vertices on the curve and t is a curve parameter. Then the derivative of the coordinates x(t) and y(t) with respect to t correspond to the orientation of the fiber. Once a point to start tracing is supplied as an initial condition, we can acquire a parametric curve for the fiber by numerical integration of the orientation field.

Fiber Tracing
As an initial value problem, our fiber tracing starts from a pixel with integer coordinates, which is confidently on the center line of a fiber. To check whether the pixel is on the center line, we follow the previous method [3] using the following two criteria: The side points P left and P right are taken from the points on the line perpendicular to the curve at P, and are sampled at the left and right side of P at a location 1 pixel away from P.
The threshold ε is a parameter to control the sufficiency of the gap between the confidences of the central pixel and its neighbors. We use ε = 0.1 by default in our implementation.
Once the start pixel P 0 = (x(t 0 ), y(t 0 )) is determined, we can trace the fiber from it by solving an ODE to generate a series of points P 1 , P 2 , . . . , P i . Note that the following tracing operation is performed using the coordinate system of real numbers, whereas the position of the start pixel P 0 is represented by integer values. Assuming that the curve parameter t is an arc-length parameter, an ODE to solve in order to extract a fiber is formulated as follows: In previous studies [3,23], the first-order Euler method is used to solve this ODE, but the ODE can also be solved by arbitrary higher-order solvers, such as the midpoint method and Heun method. We select the fourth-order Runge-Kutta method, and rewrite the above ODE in a discrete form. y 4 ), In these formulas, h denotes the step size used in the numerical integration; we used h = 1.6 in our experiments.
To compensate for the positional error accumulated through tracing, we also correct the tracing position every time after the time development step, as performed in the previous method [3]. As shown in Fig. 5, the confidence values are taken at two nearby positions P − i and P + i on the line perpendicular to the current tracing direction at P i . The points P i , P − i , and P + i are depicted by yellow, green, and blue points, respectively. Let u i = (cos θ i , sin θ i ) be the orientation of the fiber at the current tracing position P i , and u ⊥ i = (sin θ i , − cos θ i ) be the direction perpendicular to u i . Then the positions of P − i and P + i are defined as P − i = P i − Δu ⊥ i and P + i = P i + Δu ⊥ i using a constant distance Δ, where Δ = 1 is used in our system. Using the Fig. 5 Illustration of the recentering operation. In the left image, recentering is operated on the coordinates with two orthogonal directions, i.e., u i along the fiber direction and u ⊥ i perpendicular to it. A plot of the tent function is shown in the center, where the pixel P i is translated along u ⊥ i to be directly below the peak of the tent function. The result of the recentering operation for a real sample is shown on the right. In this image, the recentered fiber is shown in red confidence values of these three positions, we search for a position P * i , which is depicted by the red point in Fig. 5, between P − i and P + i where the confidence is locally maximal. This operation is performed by fitting a tent function to P i , P − i , and P + i as illustrated in the center of Fig. 5. The current position P i is moved to the recentered position P * i on the top of the tent function. More formally, the position of P * i is defined by the displacement δ * along u ⊥ i as follows: The derivation of this formula is provided in Appendix A. As shown in the right of Fig. 5, this recentering process can prevent tracing from drifting away from the thin body of the fiber.
The Runge-Kutta's time development step and recentering operation are repeated until at least one of the following two conditions is met.
The first criterion checks that the confidence at each tracing positions is sufficient and higher than w high , where C 1 = 0.7 in our experiments. The second condition guarantees the orientation of fibers changes smoothly during tracing. The tracing stops if the difference in the fiber orientations defined in (4) is greater than the predefined threshold θ max . After tracing is terminated for a fiber, we start tracing the same fiber again from the same start position, but traverse the opposite direction because the start position that satisfies the criteria in Eqs. (10) and (11) can be at the middle of the fiber. In addition, to avoid restarting the tracing from already-traced pixels, we remove the traced pixels and their neighbors with confidence greater than w high from the candidates start positions. This fiber tracing process is carried out until no start positions remain.

Fiber Segment Clustering and Length Calculation
Although the termination criteria of the fiber tracing in Eqs. (16) and (17) work well to trace a single isolated fiber, they may terminate the tracing in the middle of physical fiber near to a fiber intersection. Typical termination of fiber tracing is shown in Fig. 6. Although the tracing stops properly at the end of a fiber in Fig. 6a, it stops at the intersections due to both the insufficient confidence and sudden change of tracing directions, as shown in Fig. 6b-d. Because the purpose of our study is not only to extract the fiber locations, but also to calculate fiber lengths, these disconnected fiber segments must be connected.
To connect the corresponding fiber segments at intersections, we apply agglomerative clustering to the given set of fiber fragments. An overview of the clustering process is illustrated in Fig. 7. Let us denote an endpoint of a fiber by Q j . As illustrated in Fig. 7(i), we first collect the endpoints Q jk within a distance δ max from Q j , where k denotes an index for a point in the neighbors of Q j . Then the pair of points (Q j , Q jk ) is chosen as a candidate that can be connected when Q jk satisfies the following three conditions: The first condition (a) ensures that the two endpoints Q j and Q jk have similar orientations, which discards the endpoints The second condition (b) ensures the midpoint of two endpoints is on a fiber with sufficient confidence, which avoids connecting physically disconnected fibers when their endpoints are incidentally close to each other. This condition is illustrated in Fig. 7(ii), and, for instance, a midpoint between Q j and Q j4 is excluded from candidates because the midpoint is on a position with low confidence. The last condition (c) ensures the coherency of the fibers given by connecting two endpoints. As depicted in Fig. 7(iii), θ jk represents the change in fiber direction when two endpoints are connected. There can also be multiple endpoints chosen as a candidate pair (Q j , Q jk ) even under these three conditions. Hence, we collect the candidate pairs for all the endpoints at the beginning. The candidate pairs are stored in a priority queue in ascending order of θ jk . Then, the pairs are sequentially connected if neither endpoints in a pair is connected to other endpoints. This process of agglomerative clustering is carried out until the queue is empty. After all the endpoints are processed, we calculate the length of each fiber by simply summing up the lengths of the segments P i+1 −P i . Then, the fibers shorter than 20 pixels are discarded as outliers. In our system, we apply a moving average filter to vertices on the traced fiber before calculating the lengths, which aims to smooth the displayed fibers. Note that the smoothing operation here is just an option to make the fibers more visually pleasant and has little effect to the accuracy in fiber length measurement. The detailed effects of this smoothing operation are discussed in Appendix B.

Experiments
To demonstrate the practical applicability of the proposed system, we show experimental results for several X-ray imaged fibers in this section. In our experiment, the proposed system is implemented using MATLAB and tested on a computer with 3.6 GHz Intel Xeon E5-1650 v4 CPU and 256 GB of RAM. Because the ground truth data for the fiber lengths is not typically provided for real fiber-reinforced materials, we prepare handmade samples with reference data for the lengths and shapes of fibers. We first evaluate the system qualitatively and quantitatively using the handmade samples, and then, evaluate it qualitatively using real CFRTP samples. The hyperparameters and the values used for them in our experiments are listed in Table 1.

Experiments on Handmade Samples
The handmade samples used in the experiment are shown in Fig. 8. These samples were made by fishing lines of two different diameters, 0.910 mm and 0.520 mm, and embedded in a round piece of Styrofoam with a diameter of 300 mm. The length of each fishing line was manually measured with  a, a + b, a + 2b, . . ., a + N b, where a + N b ≤ c < a + (N + 1)b Fig. 8 Handmade experimental samples composed of plastic fishing lines as fibers embedded in a Styrofoam as a matrix a tape measure and is used as a reference value. We prepared two different samples with different layouts of fishing lines, which we named as "sample A" and "sample B," and stacked them to produce another sample named "sample A/B" (see Fig. 8). These three samples are imaged by X-ray transmission imaging using an industrial X-ray CT device (Carl Zeiss Metrotom 1500). In this imaging process, we set X-ray tube voltage to 100 kV and its current to 900 μA. Each single image is captured by an integration time of 2000 ms, and 20 images are averaged to reduce digital noise. The captured images are cropped to 900 × 900 pixels without distortion. The resulting pixel size of these images was 0.22 mm. The fiber tracing result for each sample is shown in Fig. 9. In this figure, the input radiographic images are shown at the top, and their corresponding results are shown at the bottom. The traced fibers are depicted by random colors. Further, two regions in each image are highlighted by blue and red squares. These regions are enlarged and shown to the right of the image. According to this figure, it is demonstrated that our system succeeded in tracing not only detached fibers, which we can see the blue square of sample A and the red square of sample B, but also dense intersecting fibers, which we can observe in the red square for sample A and both squares for sample A/B.
For a more elaborate analysis, we also evaluated our system quantitatively using two different standards. First, using the reference fiber bodies traced manually for each fiber, we calculated the precision and recall for the results of fiber tracing. Although these metrics are typically used to check the performance of an algorithm for problems with two choices, they can also be defined to check the similarity of two shapes [13]. For this definition, we need to use a threshold τ to determine whether a point on the estimated fiber is sufficiently close to the reference fiber. Let F be the set of positions on the estimated fibers and G be that on the reference fibers. Then, we can define the precision P(F, G; τ ) and recall R(F, G; τ ) using the following equations: where 1 {·} is an indicator function that takes 1 if the predicate in the brace is met and takes 0 otherwise. The precision and recall for the result of each sample are shown in Table 2.
For this experiment, we used τ = 2.0 in the unit of pixel to calculate the precision and recall. Considering that the size of a pixel in this sample is 0.22 mm, τ = 2.0 can ensure that the gap between the estimated and reference fibers are within the gap of the thinner fishing line with a diameter of 0.520 mm. Because the precision, recall, and their mean are greater than 0.93 for all samples, we can conclude that the traced fibers obtained by our system are close enough to the real fiber locations. However, only this result does not guarantee that the fiber lengths are also accurate, as the precision and recall may be high even when many disconnected fibers are detected at locations near the reference fibers. Therefore, we next compared the reference lengths with those estimated by our system. In this comparison, we plotted the points on a chart such that its horizontal and vertical positions are the reference and estimated fiber lengths, respectively. The charts for the three handmade samples are shown in Fig 10. As in these charts, the points are plotted near to the diagonal line, which means that the estimated lengths are close to the reference lengths. Ideally, the chart for sample A/B would be a simple composition of the results of sample A and B; however, this is not the case in our result. This is because there are more intersecting fibers in sample A/B, which makes the tracing and clustering fibers more difficult. Indeed, the gray dots in the result of sample A/B are often below the diagonal line, as indicated by a dot highlighted by a red circle. As visualized in the red square in Fig. 11, the fibers indicated in blue and light blue must be the same one, although our system failed to connect them. This happened because the endpoints of these two fibers are close to an intersection and their orientations became different due to the ambiguous orientation shown in Fig. 6d. On the other hand, some fibers become longer than they actually are, as highlighted by the blue circle. In the blue square in Fig. 11, the red fiber must be two different fibers and its right part must be connected to a blue fiber below. Because these two fibers have extremely close orientations, it was difficult for our system to detect them correctly even by setting θ max to distinguish such fibers with different orientations. Thus, our system can fail to trace some of the fibers particularly at intersections, whereas most of the fibers are traced correctly, as we demonstrated by the results so far. Later in this section, we will discuss this problem in more detail.

Experiments on Real CFRTP Samples
While we demonstrated that our system works well for handmade samples in terms of qualitative and quantitative performance, as shown in the previous subsection, it is more important to evaluate the system with real fiber-reinforced composites such as CFRTPs. The fibers in CFRTP samples can hardly be imaged by ordinary X-ray absorption radiography because the physical densities of the fibers and matrix of a CFRTP are extremely similar, resulting in insufficient contrast in a radiographic image. Alternatively, we imaged a CFRTP sample using an XTLI. The XTLI is a relatively new X-ray imaging device that can obtain three kinds of images, i.e., the absorption image, refraction image, and visibility image [14,18,27]. Among these three, we used the visibility image, which obtains significantly higher contrast between the fibers and matrix. The images we used for the following experiments are shown in Fig. 12a. They are taken by an XTLI system developed by Konica Minolta, Inc. in Japan [17]. As shown in the right of this figure, the appearances in visibility images captured by the XTLI depend on the direction of the diffraction gratings placed between the Xray source and the target object [14]. Therefore, four different appearances are obtained for a single CFRTP sample, corresponding to the four different directions of the diffraction gratings which are depicted by circled arrows at the left of Fig. 12. We refer to these images as CFRTP-x with diffraction gratings at angle x. The raw images given by the XTLI are cropped and resized to a size of 900 × 900 pixels. To demonstrate the effectiveness of our system for the fiber-reinforced composites, we compare it with the previous method by Chai et al. [3] for tracing human hair strands, which is among the most similar techniques to our method, as it uses the Gabor filtering and ODE-based fiber tracing. In this experiment, we compared our results with the previous method by three different results for orientation maps, tracing results, and FLDs, namely the histograms of fiber lengths. The results are shown in Fig. 12b-d. floatbelowskip-12pt By comparing the orientation maps, we can observe that the results of the previous method is much darker than those of ours, which means the confidences of orientations computed by the previous method are lower for the CFRTP images. Interestingly, these results using the previous method look darker than their results for hair strands shown in their paper [3] as well. We consider that these results with lower confidences are due to two characteristics of fibers in this CFRTP sample. First, the fibers imaged in the CFRTP images have various widths, while the widths of hair strands are approximately uniform. Although the hair strands with uniform widths can be detected even by setting the constant values to σ and λ of the Gabor filter kernel over the entire image, they are unsuitable for the fibers in CFRTP samples. By also optimizing these two parameters in our method, the calculated confidences are improved significantly. Second, this CFRTP sample is a RaFC and consists of short fibers with random orientations, whereas hair strands are often continuous and locally parallel. In the previous method, the values of σ and λ are set to be larger than a width of the hair strands, which helps to improve the confidences in the region where the hair strands are locally parallel. Unfortunately, the randomness of the orientations in the CFRTP hinders this method to obtain high confidence values with constant σ and λ. Because of the same reason for the first characteristic, our optimization for these parameters have improved the confidences for the CFRTP sample. Overall, our method outperforms the previous method in the points where: -our method can detect ambiguous fibers with insufficient contrast, which are hardly detected by the previous method; -our method can extract fibers with different widths, although those detected by the previous methods are mostly the same; -our method is more robust to the fibers with random orientations, while the previous method often fails to obtain high confidences.
These advantages definitely help the subsequent fiber tracing step to determine accurate shapes of the individual fibers. The third column of Fig. 12 shows the comparison of fiber tracing results. The third column of Fig. 12 shows the comparison of fiber tracing results. Obviously, our method succeeded in detecting more and longer fibers than the previous method, which owes to the diffusion process for the orientation map and the higher-order ODE solver as well as the help of the better confidences in the previous step. The effects of these experiments can be seen in the result for CFRTP-0 at the regions highlighted by blue and red ellipses. At the region of the blue ellipse, there are dense fibers with  (1) and (2) in this figure correspond to those in Fig. 10. In square (1), the fibers are disconnected and a shorter length is obtained. In square (2), the red fiber is connected incorrectly to the longer fiber segment, which resulted in obtaining a longer length At the region of the red ellipse, there are also a bundle of parallel fibers that could be only partly detected by the previous method, while our method succeeded in detecting much more fibers including unobvious ones. More importantly, the green and light blue fibers detected in the region were not connected by the previous method, as it does not perform the final agglomerative clustering step. In contrast, our method with the clustering step has successfully connected these two to obtain a continuous purple fiber. Thus, three additional processes by our method, namely the diffusion of orientation maps, application of a higher-order ODE solver, and the agglomerative clustering, have all contributed to improve the tracing results. For a more practical examination, Fig. 12d shows the FLDs calculated with four images (i.e., CFRTP-0 to CFRTP-135) of the same CFRTP sample, which are taken using diffractions gratings of four different orientations. As these four charts indicate, the FLDs calculated for the four different appearances of the same CFRTP sample exhibit similar shapes. In other words, the fibers of this CFRTP sample at the different orientations follow similar length distributions. The uniformity of the fiber lengths over uniformly sampled orientations is one of the characteristics of RaFCs, which suggests that our method can provide important insights for the inspection of the properties of real RaFCs.

Discussion
Parameter Sensitivity. Although our system obtains sufficient performance on measuring fiber lengths, the results can be affected by many parameters used in our system. Among these parameters, the results are particularly sensitive to θ max and w high . These two parameters are used in Eqs. (18), (19), and (20), which control the selection of candidate pairs of endpoints in the clustering process. When these two parameters are set to make the candidate selection stricter, the resulting fibers are likely to be fragmented, but this selection can avoid incorrectly connecting fiber segments. In contrast, when they are set to relax the candidate selection, the resulting fibers are more continuous, but more fibers can be incorrectly connected to wrong ones. Because the fibers in RaFCs have various shapes and curvatures, the selection of a set of parameters that obtain globally satisfying fiber clustering is difficult. We show several examples of how the above two parameters can affect the clustering results. As can be seen in the figure, different parameter sets for (θ max , w high ) produce different clustering results. Some of the parameter sets successfully cluster fiber segments, while others fail. For future work, we would like to investigate the optimization of these parameters for automatic and adaptive fiber clustering. Overlapping Fibers. As we mentioned in the last part of Sect. 3.1, we eliminate the start points within the diameter from the traced fiber to avoid tracing the same fiber repeatedly. However, this can also negatively affect the removal of start points on the fibers closet to other traced fibers. Two typical failure cases for that are shown in Fig. 14. As shown in Fig. 14a, two fibers marked in blue and light blue, are parallel and very close to each other. As a result, the start pixels on the light blue fiber are unexpectedly removed after the blue one is traced. This procedure to delete the start pixels also causes another problem when the part of two fibers are overlapped together, e.g., the light green and light blue fibers shown in Fig. 14b. In such cases, only one of the overlapping fibers is traced, and the overlapping parts of the other fibers Our method fails to distinguish a two close and parallel fibers, such as the blue and light blue fibers highlighted by the blue circle. It is also difficult to separate b two overlapped fibers, such as light green and light blue fibers highlighted by the red circle are lost. Unfortunately, it is inherently difficult to separate two overlapping and very close fibers only using 2D X-ray radiographic images, and we need to exploit some supporting information, e.g., that obtained by 3D X-ray CT.

Conclusion
In this paper, we introduced a computer system for measuring the lengths of individual fibers in fiber-reinforced composites, which provides helpful information for inspecting the material properties of the composites. To measure the fiber lengths, our system performs a series of processes, namely, calculating the fiber orientation and its confidence at each pixel, and extracting individual fibers by tracing and clustering. By extending the previous method [3] for tracing human hair strands, we optimized all parameters of Gabor filter kernel and succeeded in robustly calculating confident orientations of fibers for those imaged ambiguously by X-ray radiography, whereas the previous method could only extract prominent fibers. This improved calculation of fiber orientations and their confidences will aid in the fiber tracing and clustering processes. By applying fiber clustering, we have succeeded in obtaining more fibers than the previous method, although their connectivity is often insufficient particularly at the intersections and junctions of fibers. Hence, we extract individual fibers by two consecutive processes of fiber tracing and fiber segment clustering, rather than extracting the fibers only by tracing them, as was done in the previous study. As demonstrated by the experimental results shown throughout this paper, our method works significantly well due to these improvements, particularly for RaFCs that include dense, short, and frequently intersecting fibers with random orientations.
For future work, we would like to investigate how well our system works for other fiber-reinforced materials such as glass fiber reinforced plastics, glass fiber reinforced concrete, and Aramid fiber reinforced polymers, although we could not test such materials, owing to the difficult access to appropriate imaging devices for them. We would also like to explore automatic parameter adjustment in fiber tracing, which is inherently a complicated combinatorial optimiza-tion problem that is difficult to solve numerically. This would be realized if the system can score the tracing results with numerical values automatically, although we currently assess the results visually by our eyes to tune up the hyperparameters. To this end, we are interested in investigating a new cost function to evaluate the tracing accuracy.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Recenering with Tent Function
Our recentering process using a tent function is equivalent to that introduced by Chai et al. [3], but we obtain a single formula in (15). Herein, we introduce the derivation of Fig. 15 Effect of smoothing after fiber tracing. As shown in this figure, smoothing only improves the visual shape of the intersecting parts. However, it does not significantly affect the fiber lengths and distribution w/o smoothing w/ smoothing this formula. The recentering process is performed on the 2D coordinate system of a confidence axis and an axis along u ⊥ i whose origin is at P i . As we illustrated at the center of Fig. 5, the confidence w(δ) at the three points w(−Δ) = w(P − i ), w(0) = w(P i ), and w(Δ) = w(P + i ) are plotted on this coordinate system as blue, green, and yellow circles, respectively. When w(0) is larger than both w(Δ) and w(−Δ), the tent function ∧(δ) is formed by two lines: l 1 and l 2 . Here, we suppose that the central axis of the fiber is located on the top of this tent function. We therefore translate P i to P * i , whose horizontal displacement along u ⊥ i coincides with the top of the tent function.
The top position of the tent function can be obtained as an intersection of two diagonal lines l 1 and l 2 . There are two possible configurations of the three points depending on the magnitude of difference of confidences w(0) − w(+Δ) > 0 and w(0) − (−Δ) > 0. When the left difference is greater than the right one, namely when (w(0)) − w(−Δ)) > (w(0) − w(+Δ)) suffices, l 1 passes over two points w(−Δ)) and (w(0)), while l 2 passes over the third point w(+Δ)). Considering that the slopes of l 1 and l 2 have the same magnitudes but opposite signs, the equations for these two lines are obtained as follows: Similarly, when the right difference is larger than the left one, the displacement becomes The only difference between (A.4) and (A.5) is in their denominators. Because the denominator is the greater one of the left and right differences, the equations can be represented in the following single form, which is also given as (15).

B Smoothing
Before calculating the fiber lengths, our system applies a moving average filter to the vertex positions on the fibers. This moving average filter calculates an average of five consecutive positions on the fiber and updates the position of the vertex in the middle of the sequence: P k = 1 5 (P k−2 + P k−1 + P k + P k+1 + P k+2 ) . (B.7) The partially enlarged tracing results before and after smoothing are illustrated in Fig. 15. We can see that the jagged shapes of fibers are smoothed by the filter to follow the shapes of real fibers in the sample image. After smoothing, the estimated length of the orange fiber changes from 89.2 to 88.8 mm, whereas the reference length of this fiber is 93 mm. The estimated length of yellow fiber changes from 46.8 to 46.7 mm, whereas the reference length of this fiber is 51 mm. Thus, the smoothing improves the visual reconstruction result but dose not influence the precision of the length calculation significantly.