Computational Materials Science

Fiber orientation tensors (FOT) are used as a compact form of representing the mechanically important quantity of fiber orientation in fiber reinforced composites. While they can be obtained via image processing methods from micro computed tomography scans ( μ CT), the specimen size needs to be sufficiently small for adequate resolution – especially in the case of carbon fibers. In order to avoid massive workload by scans and image evaluation when determining full-field FOT distributions for a plaque or a part, e.g., for comparison with process simulations, the possibilities of a direct interpolation of a few measured FOT at specific support points were opened in this paper. Hence, three different tensor interpolation methods were implemented and compared qualitatively with the help of visualization through tensor glyphs and quantitatively by calculating originally measured tensors at support points and evaluating the deviations. The methods compared in this work include two algebraic approaches, firstly, a Euclidean component averaging and secondly, a decomposition approach based on separate invariant and quaternion weighting, as well as an artificial intelligence (AI)-based method using an artificial neural network (ANN). While the decomposition method showed the best results visually, quantitatively the component averaging method and the neural network behaved better (that is for the type of quantitative error assessment used in this paper) with mean absolute errors of 0.105 and 0.114 when calculating previously measured tensors and comparing the components. With each method providing different advantages, the use for further application as well as necessary improvement is discussed. The authors would like to highlight the novelty of the methods being used with small and CT-based tensor datasets.


A B S T R A C T
Fiber orientation tensors (FOT) are used as a compact form of representing the mechanically important quantity of fiber orientation in fiber reinforced composites. While they can be obtained via image processing methods from micro computed tomography scans (μCT), the specimen size needs to be sufficiently small for adequate resolution -especially in the case of carbon fibers. In order to avoid massive workload by scans and image evaluation when determining full-field FOT distributions for a plaque or a part, e.g., for comparison with process simulations, the possibilities of a direct interpolation of a few measured FOT at specific support points were opened in this paper. Hence, three different tensor interpolation methods were implemented and compared qualitatively with the help of visualization through tensor glyphs and quantitatively by calculating originally measured tensors at support points and evaluating the deviations. The methods compared in this work include two algebraic approaches, firstly, a Euclidean component averaging and secondly, a decomposition approach based on separate invariant and quaternion weighting, as well as an artificial intelligence (AI)-based method using an artificial neural network (ANN). While the decomposition method showed the best results visually, quantitatively the component averaging method and the neural network behaved better (that is for the type of quantitative error assessment used in this paper) with mean absolute errors of 0.105 and 0.114 when calculating previously measured tensors and comparing the components. With each method providing different advantages, the use for further application as well as necessary improvement is discussed. The authors would like to highlight the novelty of the methods being used with small and CT-based tensor datasets.

Introduction
Fiber reinforced polymers (FRP) have gained increasing relevance as engineering material in multiple sectors like automotive and aeronautic industry due to their good processability and outstanding specific mechanical properties [1]. The use of thermoplastics instead of thermosets as matrix material has become particularly interesting recently due to better recyclability [2] and faster processing [3]. Carbon fibers as reinforcement material is the more expensive but also much stiffer and stronger option compared to glass fibers [4]. However, a lot of image processing and computational methods are better developed for part, which is equipped with quantitative information about material, process and possible treatment constitutes an important idea in modern continuous, data-intensive monitoring and early detection of material behavior and failure. As an example, virtual process chain approaches (also called CAE chains) [6] can be mentioned that deal with data transferal across different scales between virtual chain and physical part. In order to quantify the convoluted microstructure of short and long fiber reinforced polymers, various properties have been developed over time. Influence on physical behavior is attributed mostly to the amount of fibers in a specific area or volume compared to the amount of matrix material (fiber volume content, FVC) [7], the length of the fibers (fiber length distributions, FLD) [8,9] and the orientation of the fibers in the material (fiber orientation distributions, FOD). Fiber orientation is generally dependent on the position in the part/plate in the long fiber reinforced thermoplastic (LFT) investigated in this work. The material flow in the underlying compression molding process influences the fiber alignment. Fiber orientation therefore has an enormous influence on the anisotropy of the mechanical properties. Although FOD are classically described by scalar distribution functions, the compact representation of fiber orientation tensors (FOT) has prevailed in process or structural simulations [10,11]. Bauer et al. described the variety of FOT in detail [12]. FOT can be determined both for specimen from finished plaques as well as for describing the initial fiber orientation state in the plastificate [13]. In certain material combination cases, quantities other than fiber orientation tensors also provide a suitable description of the orientation, such as in the case of glass fiber-reinforced sheet molding compounds (SMC). This material has a typical bundle-like structure, so that a description of the mesoscopic homogenized fiber bundles as integration in simulations makes sense [14].
From μCT images, FOT can be obtained using gradient-based methods, such as the structure tensor [15]. While FVC and FLD can also be determined experimentally by either pyrolyzing the matrix material or dissolving it chemically [16][17][18], another method (outside of calculations with homogenization models) for fiber orientation tensor identification is difficult. Furthermore, the dilemma of sufficient resolution and significant sample size arises: Especially in the case of carbon fibers which typically show diameters of around 5-7 μm (and are hence significantly thinner than glass fibers), a resolution of three to four times the diameter as edge length per voxel (which is necessary to resolve single fibers) amounts to extremely small volumes. In addition to the physical limitation of CT imaging, the scale bridging problem of microstructured materials occurs: a volume too small as specimen would result in resolving orientation too locally and these small volumes often do not represent the microstructure sufficiently well, whereas (if physically possible) a tensor from one scan of a volume too big generalizes too much compared to a meaningful averaging of many local tensors. The size of a representative volume element (RVE) is difficult to find and evaluate in the case of LFT but has been considered before in the form of statistical volume elements (SVE) [19][20][21]. Assuming small sample sizes, the determination of FOT across an entire plaque (that would allow, e.g., for the comparison with resulting FOT from process simulations or as input variable for structural simulations) would require the scan and consecutive image evaluation of a high amount of samples. However, since full-field information of FOT is highly interesting for analysis and improvement of the preceding process, an alternative to costly scan generation and evaluation was sought. Accordingly, this work is a study of the extent to which a few (in this case nine) measured tensors provide sufficient information to represent the overall plate progression. It is an interpolation question, which, however, gets its difficulty in the non-scalarity of the supporting point information. Moving away from purely mathematical or academic cases of two tensors or the like, this work explicitly aims to evaluate a benchmark of different code-based methods in the context of this engineering application. Therefore, the example used is explicitly a compilation of very dissimilar tensors.

State of the art
The fiber orientation state can be described by a probability distribution function ( ) providing the probability of finding a fiber in direction of a unit vector . The unit sphere S can be generated by the set of all possible directions of , so that The orientation distribution function satisfies three main physical conditions. Firstly, ( ) must be periodic as it is impossible to distinguish between end and start point of a fiber, hence ( ) is a symmetric function: In addition, it is normalized, i.e., the integral of ( ) over the surface of the unit sphere equals one: ( ) further fulfills the continuity condition regarding the change in with time for changing fiber orientation, which shall not be elaborated here.
Advani and Tucker [10] then simplified the description of fiber orientation states with the introduction of fiber orientation tensors which are formed by the dyadic product of the vector multiplied with the distribution function and a subsequent integration over all possible directions.
Due to the periodic nature of , all odd-order integrals amount to zero while theoretically all even-order tensors are possible; with the second-and fourth-order orientation tensors being commonly used: FOT are completely symmetric and their trace equals one. The distribution function can be recovered from the orientation tensors when they are displayed in a different form (deviatoric versions). An orientation tensor is equivalent to truncating the series by which the distribution function is given by. The distribution function is represented more accurately by higher-order tensors, yet Advani and Tucker [10] argue that the representation of second-and fourth-order tensors is usually sufficient for most applications.
The problem of tensor interpolation can be theoretically avoided for fiber orientation tensors, compared to, for example, stress and strain tensors, by switching to the scalar distribution function. Subdivision into discrete directions and subsequent interpolation would equal a Euclidean interpolation of the scalar-valued function̄. In fact, this corresponds exactly to averaging the tensor components (with weights depending on distance). This can be explained by the fact that an integral is a linear mapping. However, this standard method led to a kind of ''artificial'' isotropy, which is shown by a change of shape in the manner of a rounding in the representation form of tensor glyphs. This does not necessarily seem to represent a useful averaging, as this issue of so-called tensor swelling arose both in the field of medical technology when interpolating diffusion tensors from MRI images [22][23][24] and in the case of FOT mapping in process simulation applications [25]. This swelling effect is due to non-monotonic interpolation of the tensor determinant and the Euclidean method does not preserve the positive definiteness. It was part of scientific discussions whether a more isotropic tensor as an interpolation between two anisotropic tensors that point in different directions still constitutes a reasonable behavior of fibers. It would imply that in a region of changing flow direction first some fibers turn and others stay in the original direction until most are turned in the end. The opposite idea that most turn first a bit and later completely into the new direction would rather be indicative of tensors in the center to be not significantly less anisotropic than the two next to it. As a reference solution, the first method for tensor interpolation implemented and used in this paper is this Euclidean interpolation, which will be called component averaging (CA) from now on.
The complex Riemannian interpolation is another ''global'' interpolation method [26]. However, if more than two input arguments are used, the underlying computations can only be solved implicitly, which requires an iterative and therefore computationally expensive calculation. Since this paper explicitly seeks an application-oriented method that is as fast and simple as possible while maintaining the highest possible quality, Riemannian interpolation will not be discussed further. Another logarithmic, but explicitly solvable approach was introduced by Arsigny et al. [22]: the Log-Euclidean tensor interpolation method. This method was already considered for FOT by Krauß and Kärger [25].
A completely different approach are so-called decomposition methods. These methods make use of the fact that symmetric positive definite (SPD) tensors can be decomposed into eigenvalues and eigenvectors in spectral decomposition (cf. Eq. (13)). In terms of tensor glyphs, the eigenvalues are responsible for the shape, while the eigenvectors are responsible for the orientation of the tensor in space. Thus, shape and orientation can be weighted (according to various possible distance measures) and interpolated separately and then recomposed into a tensor. This can be done directly via the eigenvalues and eigenvectors, or via detours with the help of other invariants and, for example, quaternions. This method allows in particular the resolution of the swelling effect [23,25]. It should be mentioned here that there is not only one decomposition method, but this must be understood as a kind of umbrella term, which can be executed very differently in the individual steps. The only previous use of this concept for FOT by Krauß and Kärger [25] differs from the one implemented in this work, for example. The exact concepts used for the decomposition method will be explained in detail in the Methods chapter (cf. chapter 4) and have been previously summarized by the authors in [27], but basic considerations about the decomposition method will be elaborated in the following paragraph.
As far as the shape is concerned, a direct interpolation of the eigenvalues would be conceivable. There are also approaches which handle it this way [28]. But Ennis et al. developed the concept of orthogonal invariants [24], which seem to perform very well for physical problems. Each set of invariants decomposes the tensor shape with an orthogonal basis so that the derivatives of these invariants with respect to the (fiber orientation) tensor behave as follows: Ennis et al. established the so-called K-and R-invariants, which are also used in this paper. They argue that while eigenvalues form a set of orthogonal invariants as well, they have the disadvantage of not isolating essential attributes of tensor shape like size and anisotropy which the sets of K-and R-invariants provide instead.
A direct interpolation of the eigenvalues has still been performed as a test but indeed seemed to distort the results and will not be considered in more detail in this paper.
When it comes to the interpolation of the orientation of the tensors, there are multiple works which focus on the interpolation of rotations in 3D, i.e., elements in the 3D rotation group, which is also called (3), as a separate mathematical problem independent of tensors or the shape of the same [29][30][31][32]. However, many of these papers focus on the smooth rotation between two or few different orientation states. These methods do not necessarily perform just as well for multiple and spreaded orientation tensors. Typically, the orientation of a general basis in 3D linear algebra is described by Euler angles with respect to a fixed coordinate system. Generally, it must be considered that an orientation can be described by 24 different coordinate systems. This ambiguity is counteracted with the help of conventions, this includes the determination to use a right-handed system as well as the sorting of the eigenvalues according to magnitude. However, after the conventions still four possible coordinate systems remain to describe an orientation. Depending on the choice it can be influenced whether between two orientations, which, e.g., lie only 20 • apart, the interpolated tensor rotates by 10 • , which corresponds to the -at least at first -obviously correct option, or by 170 • , which would correspond to the mirrored coordinate system. For two tensors between which one wants to interpolate, it therefore makes sense to implement a query and restrict the size of the angle. However, for a set of measured tensors T at multiple supporting points, where the ones further away from the one that is currently to be calculated have less weight but are still included in the calculation of this interpolated tensor, this becomes less obvious. This aspect will be taken up again in the discussion.
A much-investigated method is the orientation interpolation via quaternion. Quaternions as described by Hamilton extend the complex number system and are usually represented in the following form: The result is a four-dimensional number system (mathematically: a vector space) with a real part consisting of one real component and an imaginary part consisting of three components, which is also called the vector part. Multiplication of quaternions is noncommutative. Quaternions allow in many cases a computationally elegant description of three-dimensional Euclidean space, especially in the context of rotations. By using quaternions instead of Euler angles the problem of Gimbal lock can be avoided and they are simpler to compose. Compared to the rotation matrices, quaternions are more compact, efficient, and numerically stable. Regardless, not all ambiguities can be circumvented. The unit quaternions and − represent the same rotation. This means that there is a 2:1 homomorphism from quaternions of unit norm to the 3D rotation group (3). In other words, (3) is double-covered by quaternions. This sign ambiguity has to be paid attention to when computing a quaternion from the rotation matrix.
Alternative fiber orientation interpolation methods were developed, e.g., by Köbler et al. [33]. They developed a mechanical interpolation method in a surrogate model, i.e., they first calculated the material response for discrete fiber orientations and then used linear interpolation on the fiber orientation triangle (a material model as a function of the orientation triangle). While this yields good results, it is not related to the interpolation of orientations on component level and is therefore of no interest in the application thought of in this work.
The use of artificial intelligence (AI) for fiber orientation tensor interpolation has been explored by Sabiston et al. [34]. The authors used a large number of FOT obtained -as in this work -from μCT images of multiple plates of the same process. This represented their ground truth, which they used to train the artificial neural network (ANN), which was subsequently able to predict the tensor components for plates of this process with less deviation than the variability was between neighboring microstructure units. Since this represents an entirely new, non-physical, nor classically linear-algebraic way, a neural network is used as the third interpolation method in this work. However, in contrast to Sabiston's work, the ANN is trained using only the nine measured tensors considered for all methods for comparison. Then, the remaining 160 are predicted using the trained network. Even though the nine tensors give six values of information each, the use of AI with such a small amount of input data is rare, but it is intended to assess whether this can still produce reasonably useful results, or whether a useful result can be expected with a small additional number of given tensors.
Consequently, this study focuses on determining a full-field distribution of fiber orientation tensors across an entire carbon long fiber reinforced polyamide 6 plaque by interpolating orientation tensors determined from small samples at specific positions. Therefore, nine samples of 10 mm × 10 mm × 4 mm were scanned in the CT and the orientation tensor was determined for each scanned volume via structure tensor approach implemented by Pinter et al. [35]. Following, the nine computed tensors of second order were interpolated via three different methods to predict 160 additional fiber orientation tensors inbetween. The results were compared visually based on their physical meaningfulness. Subsequently, they were compared quantitatively by leaving one of the measured tensors out of the input data provided to each method and computing the originally measured tensor instead. The difference between the computed and the measured tensor is evaluated through the deviation of the absolute component values and via the Frobenius norm. The implemented interpolation methods include Euclidean interpolation (component averaging), a decomposition-based method interpolating the shape of the tensors with the help of orthogonal invariants and the orientation of the tensors using quaternions and a method based on artificial intelligence (AI) with an artificial neural network (ANN). While most of these methods have either been developed from a mathematical or theoretical point of view or -if a physical use-case existed -worked with a high amount of input tensors (in case of the ANN), the application of these methods to sparse but measured CT-based orientation data has not yet been explored to the authors' knowledge.

Notation
Symbolic tensor notation is preferred in this work. Scalar values are denoted by standard Latin and Greek letters, e.g., , . Tensors of first order are represented by bold lowercase letters, e.g., , and bold uppercase letters are used for tensors of second order such as , . Fourth-order tensors are denoted by double-struck letters like C, S.
Sets, i.e., collections of quantities, are denoted by calligraphic symbols, e.g., A and are constructed by curly braces. In them, the elements typically are given explicitly or expressed by conditions to be fulfilled by each element contained in the set. The special orthogonal group (3) represents all 3D rotations. Four-dimensional quaternions are represented by an arrow-head above the Latin letter, such as in ⃖ ⃗.
The terms () and () are the trace and determinant operators respectively, ‖ ‖ represents the Frobenius norm of the tensor defined The rotation of a tensor is denoted by the Rayleigh Product ⋆.

Scan acquisition and determination of fiber orientation tensors
The authors used a YXLON-CT precision μCT system with a flat panel PerkinElmer Y.XRD1620 detector with 2048 pixel × 2048 pixel for the acquisition of the microstructural images. The parameters of the performed μCT scans are listed in Table 1. The resulting volumetric images are reconstructed applying the Feldkamp cone-beam algorithm [36].
Subsequently, the reconstructed scans were processed in VG Studio Max 3.4.2 and cropped into regions of interest (ROI). If necessary, brightness and contrast were then adjusted in the ImageJ (FIJI) software. In addition, the individual gray value threshold was determined for each scan, which was required for the use of the method by Pinter et al. [35] that determined the FOT from the scan data. Pinter et al. implemented their code in C++ with the help of the ITK library. The algorithm makes use of the structure tensor (cf. Eq. (8)) as it performed best out of three different implemented methods: with being a discrete function describing the gray-value intensity of the reconstructed image. As the structure tensor is based on the first numerical derivative, it is known to be a robust approach. In the used algorithm, the structure tensor calculation is combined with a Gaussian blur of a width of . In this work = 0.2 voxels was chosen as it showed meaningful results across all scans. Another parameter that had to be adjusted before running the algorithm was the mask size. The obtained tensors were averaged by another Gaussian blur representing the image region that is taken into account to determine an orientation at a certain point. It has to be larger than the first one that was used for the derivative and was chosen to two voxels in this work. The algorithm then calculates the smallest eigenvalue and its corresponding eigenvector of the resulting tensor as fiber orientation. The FOT calculated with this algorithm constitute the foundation or ground truth that are fed into the different interpolation methods. They are henceforth called ''measured values'' implicitly including thathowever -these FOT are subject to a certain error as well.

General
In general, all interpolation methods were mainly implemented in Python 3.8.
Symmetric positive definite (SPD) tensors can be visualized as superquadric tensor glyphs [37,38]. This method will be used in this work as it constitutes a descriptive and interpretable way of assessing the success of the different implemented interpolation methods. The authors implemented the visualization in Matlab R2020b with the help of the ''plotDTI'' function of the fanDTasia ToolBox by Barmpoutis et al. [39].
The overall idea of all three implemented interpolation methods is to get FOT values for the 160 positions in the plaque that are missing, from calculations with the measured nine FOT at the given grid As a weight function, multiple options are conceivable with the possibly simplest being Shepard's inverse distance weighting (IDW) as an explicit approach with denoting an arbitrary point that shall be interpolated, being a known interpolating point and d is the given distance from the known point to the unknown point . is a positive real number, called the ''power parameter''. Weight decreases as distance increases from the interpolated points. Greater values of assign greater influence to values closest to the interpolated point, which results in nearly constant interpolated values for large values of .

Orientation distribution function (ODF) and component averaging method
Recalling the definition of second-order orientation tensors as described by Advani and Tucker [10],Kanatani [11], with S being the unit sphere and d the surface element on it, as well as being the unit vector for the direction of the fibers, it appears that is linear in ( ). Assuming the surface can be divided into two equally sized areas S 1 and S 2 with two distribution functions 1 ( ) and 2 ( ) and ( ) = 1 2 ( 1 ( ) + 2 ( )) holds, this means that = 1 2 ( 1 + 2 ) is exact, as integration is a linear functional and as an integral domain can always be divided into sub-intervals. This further implies that a direct averaging of the orientation function is equivalent to an averaging of the components of the orientation tensors.
Thus, the algorithm multiplies the components of each measured FOT by a weight that depends on its distance from the tensor being calculated.
As mentioned before, Shepard's inverse distance weighting method is used as weight function in all three methods with = 2: Compared to Eq. (10), Eq. (12) features a necessary normalization factor.

Decomposition-based method
The method, which uses spectral decomposition of tensors, is shown schematically in Fig. 1 and is described in the following.
For the chosen decomposition approach, the shape and orientation of the tensors are to be interpolated separately. Therefore, the wellknown spectral decomposition resulting from the Eigenvalue problem is used: denotes the tensor containing the Eigenvalues on the principal diagonal and is defined as the orthogonal rotation matrix consisting of the normalized Eigenvectors. Orientation The rotation matrix R can be interpreted as a rotation around a rotation axis and can therefore be transformed into a quaternion as described in the state of art: with rotation axis: = ( , , ) and rotation angle . The quaternion is calculated from the given rotation matrix via: This is followed by the actual interpolation: = ∑ with weights: ∑ = 1 and the retransformation of in : Projector interpolation An alternative approach for orientation interpolation, circumventing the sign-ambiguity, is based on the bi-orthogonal projector decomposition with Gahm and Ennis [40] introduced this method as Dyadic interpolation, since each projector can be expressed as dyadic product of its underlying eigenvector . Subsequently, the projectors are interpolated separately and composed into a non-orthogonal deformation tensor . Finally, polar or singular value decomposition yield the interpolated orientation tensor . Details on the implementation can be found in [40]. It is noteworthy, however, that this approach will fail, if the input projectors are antipodal.
Shape For the interpolation of the shape, three linear independent invariants are formed of each tensor and interpolated separately. Of the already mentioned K-and R-invariants [24] 1 , 2 , and 3 will be used (comparable to [23]): Even though 1 and 2 are not orthogonal (cf. appendix of [24]), the use of 1 can be justified by ensuring that the trace of the orientation tensor is one. It is not necessarily essential to have orthogonal invariants for this specific case of application.
The invariants are then interpolated individually: The weights stay the same: ∑ = 1. From the interpolated invariants, the following formula was used to calculate the associated eigenvalues (cf. [23]): For = 1, 2, 3 holds: With these eigenvalues, can then be created again.

Artificial neural network (ANN) method
The artificial neural network used in this study is based on the idea and implementation of Sabiston et al. [34]. Just like for the other two methods, the goal of the neural network is to determine a FOT for each specified point of the 160 missing positions within the plaque. For the ANN, the nine x, (and z) coordinates of the given FOT were normalized (divided by 14 since there are 13 rows and columns of FOT). Thesẽ,̃, and̃represent the input data. The output data for the ANN are the respective components of the nine given orientation tensors of second order. Since these components are already between −1 and 1, this data does not need to be normalized. The coordinates were read in as one .csv file as input and the components separated by 11 , 33 , 12 , 13 and 23 as five separate .csv files as output. There are only five independent components instead of the usual six independent ones for symmetric tensors, since orientation tensors are subject to an additional condition that the trace of the tensor must add up to 1: Therefore, only 11 and 33 were fed into the network as output parameters and 22 was determined via 22 = 1 − 11 − 33 . The choice was made explicitly, according to the findings of Sabiston et al. [34], to use only one in-plane coordinate and the through-thickness coordinate ( 22 and 33 would have worked analogously as well) in order to reduce the error and to satisfy Eq. (19). This is due to the two in-plane coordinates being significantly larger than the throughthickness coordinate, which in turn meant that 11 and 22 alone could get above 1 if they were both predicted.
The ANN consists of an input layer, where the normalized input coordinates of the = 9 different points are given and two hidden layers with = 48 neurons. The output of the first hidden layer is the input for the second hidden layer. In the output layer the five independent tensor components 11 , 33 , 12 , 13 and 23 are predicted for the given = 9 points. The structure of the ANN can be seen in Fig. 2.
The optimizer used is the stochastic gradient descent (SGD) which is an iterative method for optimizing an objective function ( ) with suitable smoothness properties. Thereby, after choosing an initial vector of parameters and a learning rate , two steps are repeated until an approximate minimum is obtained: The samples in the training set are randomly shuffled and ∶= − ∇ ( ) is set for = 1, 2, … , (with being the summand function typically associated with the th observation in the data set (used for training)). The loss function chosen is the Mean Absolute Error (MAE), which is defined as follows: with being the number of samples, ( ) being the value of the orientation tensor component at the sample location and̂( ) being the predicted value of the orientation tensor component at that sample location. MAE was preferred as error metric over percentage error since many values (especially the off-axis and 33 components) are close to zero. Therefore percentage errors tend to become quite large. Additionally, outliers seem to be filtered out better by using MAE than by a quadratic error metric like root mean squared error which is more likely to result in overfitting and being biased towards outliers respectively. Additionally, a soft sign activation function is used in the model as it is able to calculate negative numbers and behaves differently in terms of saturation (compared to, e.g., the hyperbolic tangent) because of its smoother asymptotes (polynomial instead of exponential) [41]. However, this of course impacts the amount of epochs required for training as it does not saturate as quickly. The soft sign activation function is given as where is the input to the function and is the output of the function. Furthermore, a bias was placed on the loss weights (w) of the outputs of the ANN in order to give more weight to the in-plane orientations. The biases are 0.4 for 11 , and 0.15 for all four other components, adding up to 1.
The classical data validation split of 25% is used in the study. The high amount of 100.000 epochs, i.e., times the neural network iteratively trains the weights for each neuron to optimize the outputs from the given input steps, was chosen. While increasing the epochs normally reduces the error, it can also evoke overtraining and leads to longer calculation times. The authors will elaborate on that in the discussion.
All chosen parameters of the ANN are summarized in Table 2. Once the model is trained, a .csv-file with all 160 normalized coordinates -apart from the nine the network was trained with -is given to the trained ANN, to predict the components of the missing FOT.

Material
In this study, carbon fiber reinforced polyamide 6 (PA6) is investigated. This material is manufactured in the so-called ''long-fiber thermoplastic direct process'' introduced by Krause et al. [42]. The schematic process can be seen in Fig. 3.
The LFT-D process is characterized by its inline compounding of the matrix polymer (in this work PA6) in a regular extruder at about 230 • C. The polymer melt is then fed through a film-die into the second extruder which is a twin-screw-device and which is used as mixing unit. The (through a specific nozzle) ejected plastificate is subsequently manually inserted in a Dieffenbacher compression molding press which prevents a reheating of semi-finished products that is necessary in other known LFT production processes. It is then pressed and cooled in the tool which has dimensions of 400 mm × 400 mm and a temperature of about 80 • C −90 • C. The variation of the plastificates in height and size as well as the at times uneven temperature in the tool, which was measured by the authors mid-process, might be important contributors to variations in the flow behavior of the melt and therefore the fiber orientation behavior as a function of position. This aspect is elaborated in the discussion.   Out of one of these manufactured plaques nine rectangular samples of 10 mm × 10 mm × 4 mm were cut via waterjet cutting (cf. Fig. 4).
The right-hand coordinate system in the lower left corner of the plate is also valid for all following images with tensor glyphs. It is shifted slightly to the upper right for better visibility, but should actually have its origin at the outer corner of the plate. The dimensions of the small samples resulted from the consideration of two main aspects: How small do the samples have to be in order to allow a sufficiently fine resolution of the fibers and how large do they have to be in order not to cut off the long fibers too much and thus possibly distort the result of the detectable orientation tensors. This question plays somewhat into the aforementioned conflict point of RVE in the area of fiberreinforced or generally heterogeneously microstructured materials. The first question was clarified by a series of experiments: several scans with different resolutions and therefore sample dimensions were made and generally compared in their image quality. While single fiber detection (which is required for example for image recognition of fiber length distributions) requires an amount of 3-4 voxels over the fiber thickness and thus in the case of carbon fibers approximately a resolution between 1-2 μm/voxel, the fiber orientation algorithm works on voxel data without any connectivity. Hence, it is not necessary to separate fibers and matrix in a pre-processing step and a resolution of one voxel over the thickness is reasonable. As described above, a resolution of 8.57 μm/voxel was chosen in this work. The second question was answered by experimentally determining the fiber length distribution in the plaque. A peak at about 500 μm and maximum lengths of about  15 mm resulted in a weight-averaged length of 1.66 mm and a valueaveraged length of 0.4 mm. The sample area of 10 × 10 mm 2 could thus be well justified, since the fewest fibers would be cut off.

Fiber orientation tensors
First, μCT scans of the nine samples were made as described above. An example 2D image from the center of such a scan (in this case of the sample UR) is shown in Fig. 5.
Then, one FOT per sample was determined from these scans using the structure tensor method described above. The result of the structure tensor based evaluation of the μCT scans resulted in the following FOT for the set T = { , , , The nine tensors, which are displayed in Table 3, are therefore the basis of all further calculations. They are only one possible example of FOD of a single plate, but they represent a rather extreme case of fiber orientation in such a plate, which is why the different behavior of the individual interpolation methods can be shown relatively well. In Fig. 6, these tensors can be seen represented as tensor glyphs at their respective positions in the plate.
The tensors in the areas in between (marked with a question mark) are to be determined from the nine given by means of the different methods.
It should be mentioned that the measured tensors were determined over the complete sample volume and were thus also averaged over the thickness. However, especially over the thickness, the tensor components still change strongly even in such a small volume. In Fig. 7, where the principal components of all fiber orientations tensors are plotted over the thickness, one can see strongly heterogeneous courses.
In particular, the change of orientation in the edge layers can be detected in all tensors. In the left and upper tensors, the 22  This edge phenomenon, known in compression molding as the ''shell-core effect'', results from the rapid freezing of the hot fiberpolymer mass (about 230 • C) at the ''cold'' (about 80 • C −90 • C) mold wall, while the central region longer remains fluid. This central layer flows forward in a so-called swelling flow from the center into the not yet filled mold area. Fig. 8 explains, why at the position of , the 33 component does not increase at the edges of the height of the plate. The tensor was determined from a specimen located in the area where the plastificate is inserted (cf. Fig. 4). In the μCT scan of the plastificate in Fig. 8, the upper and lower region show pronounced orientation in direction (hence ''11-direction''), which is transferred onto the plate as the mass freezes directly with the closing mold in these areas (the plastificate is much thicker than the final plaque). Therefore, for , the 11 component increases at the borders (cf. Fig. 7. On the contrary, one can see the increase in the 33 component for and in Fig. 7, which results from the mass flowing upwards and downwards at the edges of the height (Fig. 8). For the position, this occurs especially pronounced as the mold ends and the movement of the mass in both the x-and the -direction stops, resulting in increased flow towards both positive and negative -direction from the center of the mold.

Component averaging method
The set of measured orientation tensors T via CT scan and subsequent calculation via structure tensor is represented by the blue tensor J. Blarr et al.   Fig. 7. A real μCT image of one of the plastificates used in the compression molding process can be seen at the left. glyphs in Fig. 9, the set of interpolated tensors T = { ∀ ∈ 1, … , 13 ∩ ∈ 1, … , 13} by the orange tensor glyphs.
The origin of the global coordinate system is located in the lower left corner of the plaque. The original LFT charge covered almost the entire left side of the 400 × 400 mm 2 mold with a width of about = 90 mm (to the right), a length of about = 350 mm (up) and a height of about z = 60 mm. Thus, when the press closes, one would expect a quasi 1D flow to the right.
However, in Fig. 9 a clear progression can be seen in the measured (blue) fiber orientation.
After a clear preferred direction in the left region resulting from the plastificate, i.e., from the last extrusion step in the LFT-D process, a turn to a rather dominant orientation in the -axis seems to happen towards the middle of the plaque (apart from the top region). At the right side another turn to an again more -direction-dominant orientation happens (apart from the lowest tensor glyph, which in general seems to be more isotropic than the other tensors). Considering these measured tensors, the interpolated tensors should follow some kind of curve. In fact, the interpolation does not seem to cover this behavior smoothly but instead rather accomplished the orientation changes through rounding the tensors. Following the literature, this behavior was expected (cf. swelling effect in the State of Art (chapter 2) and can be confirmed.
In order to be able to approach quantitative error analyses and to better assess the interpolation behavior, one measured tensor of T was omitted in each case and instead also determined with the interpolation method. The visualization results are shown in the nine lower pictures in Fig. 10.
To obtain a quantitative error value, the Frobenius norm of the measured tensors and their respective interpolated substitutes was formed. The result of the difference between the Frobenius norm of the interpolated and the original tensor can be seen as an error map in Fig. 11. The method seems to predict the , and tensors the worst. It is difficult to judge whether the Frobenius norm is suitable as a quantitative assessment, but it will be discussed further in the discussion.
Therefore, a third possibility of error analysis is considered, namely the direct component comparison between interpolated and measured tensors. Fig. 12 shows this for the component averaging method.
The nine graphs in Fig. 12 correlate to the nine different tensor components of a 3 × 3 tensor. Each graph shows the component value of the measured tensors in blue and of the interpolated tensors in orange on the -axis for each of the nine measured tensors (depicted on the -axis, starting with ). While the differences are largest for

Quaternion-based
In Fig. 13, the results of the interpolation with the quaternion-based decomposition approach can be seen.
The before-mentioned progress of orientations can be visually traced as a clear curve. As for the interpolation method as such, the visual results are for the most part very appealing. Interpolation between the individual measured FOT is good and the transition between two adjacent tensors also appears reasonable. The anisotropy is not basically lost between two differently oriented tensors by ''rounding the tensor''. The rotation of two adjacent tensors occurs with small angles and therefore smoothly. The only exception to this can be seen at the upper right edge: The interpolated tensor 10,13 in the middle of and behaves somewhat strangely as far as the behavior of the row is concerned. Instead of closing the estimated angle of 20 • between and by a piece wise change of 10 • , the tensor 10,13 is oriented in an angle deviating by around 80 • compared to the measured ones. However, the tensor is for example also taken into account for the calculation of this tensor (just like all the other measured ones of the set T ), even if weighted less strongly than and , which favors the big rotation of the tensors in the uppermost row considering the global orientation behavior. Furthermore, the behavior in this column looks much better than could be expected if the tensor had rotated in mathematically negative (clockwise) direction around the -axis than the chosen positive (anti-clockwise) direction. When leaving measured tensors out of the ''ground truth'' and interpolating them instead, there are definite changes in the orientation course, which can be seen in Fig. 14.  Fig. 13. Visualization of interpolated (orange) and measured (blue) tensors when using the decomposition-based interpolation method described in this paper.

Fig. 14.
The graphic shows the measured (blue) and interpolated (orange) fiber orientation tensor glyphs when leaving one measured tensor out of the calculation and interpolating it instead with the decomposition method respectively.
For example, the behavior of the afore-mentioned 10,13 changes significantly when or are omitted. In general, however, it must be stated that the orientation courses react very agilely and sensibly to the changes when individual tensors are omitted when using the decomposition method.
The quantitative evaluation based on the Frobenius norm is visualized in Fig. 15. The rather poorer interpolation at the left and upper edges and the relatively good performance in the middle of the plaque (and lower right) are noticeable.
Considering the component-wise deviations displayed in Fig. 16, it is striking that they are considerably high in this specific case for this  coordinate system. However, it is noticeable that major trends between the different tensors are mostly preserved by this method (cf., e.g., 11 ).
Projector-based For the sake of completeness, projector-based interpolation has been applied to the experimental data. The according results are displayed in Fig. 17.
Analogously to the quaternion-based approach, tensor shape is interpolated monotonously. Thus, the interpolation does not yield results with larger isotropy than included in the data set. However, significant deviations can be observed regarding tensor orientation. While the rowwise ''flipping'' behavior appears to be resolved (compare first and last row of 16 and 17), a singularity right next to the domain's midpoint occurs. As the implementation of the projector-based interpolation is based upon iterative methods, it cannot be vectorized. Therefore, the authors decided to neglect this method and concentrate on the quaternion-based approach for all further investigations in this work.

ANN method
The results of the tensor field when interpolating with the neural network can be seen in Fig. 18.
It strikes that when training this network with the measured tensors, it is able to produce both very anisotropic and very isotropic tensors at the points with missing tensors, compared to the other two main methods that rather dispensed one or the other. While some areas look smooth like, e.g., the upper and right area, there is non-monotonous interpolation behavior concerning, e.g., and where the surrounding tensors are much more isotropic than the measured one and Fig. 17. Visualization of interpolated (orange) and measured (blue) tensors when using the projector-based interpolation method described by Krauß and Kärger [25]. Fig. 18. Visualization of interpolated (orange) and measured (blue) tensors when using the ANN-based interpolation method described in this paper. also quite differently oriented (looking specifically at the tensor ). When training the network multiple times with the same input data, the results look very much alike, speaking for the robustness of the method.
However, looking at the plots of the tensor fields when leaving measured tensors out of the input data (Fig. 19), some of the nonmonotonous behavior shows again (cf., e.g., the fields without and without ).
Consequently, the difference of the Frobenius norm for these two cases is significantly large and so is the error for the interpolation of the tensor, which is all depicted in Fig. 20.  Considering the component-wise errors in Fig. 21, the results are still comparatively good. The largest deviations occur especially for the three mentioned tensors before. Just as was the case for the other two methods, the biggest deviations appear mostly for the 11 and 22 components. The comparison between the three methods regarding the differences between the measured and interpolated tensor components will be taken up again in subchapter 7.2.
The quality of the neural network must also be considered with respect to the progression of an error measure over the number of epochs. As an error measure, the already mentioned MAE was used and the course over the epochs can be seen in Fig. 22. As expected with the small number of training data, the network does not behave optimally. However, the graphs show both overfitting ( 23 , 33 , 31 ) and underfitting ( 11 ) trends. Thus, it is difficult to draw a general conclusion. Overfitting rather argues for using more training data or stopping at a lower number of epochs, apart from solutions that require more specific methods which are very dependent on the model. Underfitting can be combated by different measures depending on the cause of the underfitting; either more epochs (unlikely here) or more parameters of the model can help, or a change to a completely different model.  However, more training data can also help with underfitting, which is most likely in the case considered.

Curve in fiber orientation
Although this work only shows the fiber orientation field of one plaque, scans and also tensile tests in different angles of further plates hint towards a slight drift of the mass to the top too. Uneven height or dimensions of the plastificate and uneven temperatures in the tool could explain the curve in the orientation tensor field. The authors distinctively measured slight variations in both temperature of the mold and geometry of the plastificate. Furthermore, the initial charge is inserted in the mold manually, resulting in a possibly skewed/angular position of the latter when the press closes. However, the authors assume the main reason for the drift to be the ejection of the initial charge from the extruder. As this process takes a few seconds, the ''old'', first ejected plastificate has been cooling down for longer than the part that leaves the extruder at the end. The colder part is located in the upper part of the plaque (higher value in the coordinate system given in this work). As a result, the lower part of the fiber-matrix-mass should be less viscous and better flowable, presumably causing the flow front to advance faster there and hence eventuating in the curve.

Comparison of the three methods
Looking at the deviations of the measured and calculated components for all three methods in Fig. 23, it is noticeable that the component averaging method actually performs best for this particular case.
However, the question remains whether an averaging can physically correspond to the same idea as that of an interpolation. After all, the search is not for the ''mean'' value, but for the ''middle'' of two given tensors. The already mentioned difficulties of this method concerning the shape of the interpolated orientation tensors rather pleads against this method, but the frequency of use and the simplicity of implementation and calculation are of course advantages.
The neural network is also very close to the component averaging method in terms of absolute deviations. The neural network is explicitly trained on the components of the measured tensors and tries to predict them, so this is not surprising.
The decomposition method performs worst with respect to the value deviations of the components. Table 4 sums up the mean average absolute error of the different methods.
However, it has to be stated that the retroactive computation of originally given tensors does not correspond exactly to the task of the determination of tensors between given ones. The former works with an input data set reduced by 1/9. In particular, when leaving out one of the four corner nodes, rather the extrapolation instead of the interpolation behavior of the respective method is demanded and judged as the calculated point is located outside of the support point grid. In fact, when calculating the average absolute error of the different methods without including the prediction of the corner tensors, the mean average error gets significantly smaller as can be seen in Table 5.
The deviation to the quantitative error when involving all nine tensors is the highest for the decomposition method. This could lead to the assumption that the decomposition method performs so poorly quantitatively because this error assessment is not appropriate for the   [25] as well as mathematical accuracy clearly speak for the superiority of the decomposition approach. This is supported by the fact that the decomposition method performs significantly better without the extrapolation task (omitting corner points). Additionally, the reduction of the input tensors is a change of the initial situation and only slightly resembles the actual interpolation task. So, judging the methods in different ways is necessary. While concerning the component deviations the decomposition method performed the worst out of the three methods, the visually agile interpolation behavior in the tensor field has already been mentioned. The authors conducted a study, where the measured tensors were expressed in different coordinate systems (rotated by different angles). For these new tensor components, it was remarkable how well the decomposition method interpolation adapted itself. Furthermore, the monotonous behavior concerning anisotropy is preferable compared to the other methods. In addition, the authors highly assume that the decomposition method might perform better when the tensors are located further in the center of the plate/part, since the decomposition method interpolated worse at the edges. This will be taken up again in the outlook.
The two different versions of the decomposition method differ in the use of projectors instead of invariants. The former have the advantage that they are unique for a given tensor. While this prevents the effect observed at the upper right edge of the plaque of the invariant decomposition method, it shows a singularity close to the domain's midpoint. The visual smoothness of rotation depends on whether one considers the row or the column course of rotation. The question remains whether the measured tensors of this specific plate are generally too inhomogeneous to represent the original fiber flow optimally. After all, the sample sizes of 10×10 mm are small and can also map very local effects. Again, this will be addressed in the Outlook.
Considering adaptability, the ANN provides the significant disadvantage of being trained for a specific grid of measured tensors. So the network trained on those will not perform well for a different plaque. However, this work was precisely about the quick characterization of a single plate or part. For this rather unusual application for a neural network, the possibilities to achieve meaningful results relatively fast with little data should be discussed. With run times of less than one hour for a network with two hidden layers as well as for one with three (which has also been tested by the authors), one can definitely speak of a time saving in comparison to a one and a half hour scan and additional computation time for the generation of a tensor via μCT scan. As far as the component comparison between measured and calculated tensors is concerned, the neural network behaves similarly well as the component averaging method. However, both decomposition method and component averaging are both significantly faster as they generate results near-instantly. The neural network is more suitable for a similar application: Just as in [34], the fiber orientation field can be predicted for one specific, stable process by evaluating multiple plaques via μCT scan and subsequent FOT determination and training the network with a high number of FOT. This provides for good results of a general fiber orientation of a stable process.
In general, the problem of the issue addressed in this paper lies in the combination of the mathematical implementation of orientation and rotation in three-dimensional space and the physical effects of the flow of material. In particular, the assessment of the extent to which an interpolation is visually or quantitatively convincing is also inevitably related to the question of how many mathematical bounds one may give so that the tensor interpolation becomes meaningful, but one does not artificially interpolate in such a way that it ''looks better''.
While neural networks usually ignore any implementation of physical effects (except for physics-informed neural networks (PINNs)) and learn only from previous data, the other two analytical methods also do not consider boundary conditions (such as the location of the initial charge in the die in the press), although this would be theoretically possible. On the contrary, these are also purely mathematical considerations of tensor interpolation, since a large intervention in the interpolation ought to be avoided and the main focus was on the given data. The aim was to find a relatively easy to implement and relatively well-performing method for fast interpolation of given orientation tensors of small number to avoid costly scan times.
In fact, Brannon even argues in chapter 14 of her book ''Rotation, Reflection, and Frame Changes -Orthogonal tensors in computational engineering mechanics'' [43] that methods for mixtures of rotations must be selected appropriate to the physical application. The motivation for this paper therefore lies in finding a reliable tensor interpolation method for this specific application case. For any other application than fiber orientation tensors, this benchmarking study might have to be repeated (even though the decomposition method seems to perform quite well for both diffusion and orientation tensors). In addition to the application, the amount of data or grid points also might influence the functionality of the specific method. Many original algebra-derived methods might work for two given tensors with limited tensors to be interpolated in-between, but perform way worse for many given tensors or many tensor to be calculated. Conversely, many given tensors could certainly improve the result of the methods presented here.

Conclusion
In summary, the authors applied four different tensor interpolation methods to fiber orientation tensors. The implemented methods provide for good macroscopic interpolation results for fiber orientation tensor fields interpolating measured FOT determined from microscopic Xray computed tomography scans. The authors were able to generate orientation fields of plaques in significantly less time than scans for all of these positions would have taken. A particularly interesting aspect is the use of only nine input tensors, showing that the presented approaches work for scarce data. These methods could therefore replace high amounts of resource-consuming μCT scans. Each method has its advantages and disadvantages. The observed monotonous behavior of the decomposition method is particularly satisfactory. When calculating originally measured tensors, the performance of the other two methods were more convincing. The mean absolute error when calculating measured tensors ranges from 0.105 for the component averaging over 0.114 for the neural network to 0.204 for the decomposition method using invariants. However, as this assessment includes an extrapolation task and a changed initial situation, this cannot necessarily be attributed to the quality of the interpolation procedure itself.

Outlook
There are multiple possible adjustments of the methods shown in this work. The interpolation behavior of the invariant decomposition method at the top right could be artificially prevented by restricting the possible angle. This could be implemented by comparing the results of the scalar products of the specific quaternions involved and choosing the combination providing the maximum scalar product and hence the smallest possible angle. However, the question arises whether this is not an intervention, since the influence on the interpolation of the tensors in the immediate vicinity would then be significantly higher than given in this work. The authors tried to implement as little artificial restrictions as possible. Furthermore, is not normalized in the invariant decomposition method. And if several unit quaternions are added, the result is no longer a unit quaternion. Analogously, this is true for rotation tensors. Thus, the rotation interpolation as implemented in this work also changes a bit the shape of the tensors. By a normalization a further adjustment could take place here.
However, the most important and perhaps most obvious next step according to the authors -after benchmarking different methods for this interpolation application in this work -is a sensitivity study regarding plaques, sample geometry and sample location. The performance of the methods should definitely be studied for different orientation fields, i.e., different plaques. In a small way this has already been done, but the focus of this paper was rather on the difference between the methods for a specific case. Furthermore, it would be important to take more and also less samples from a plate and evaluate which amount is too little to result in a meaningful orientation field and vice versa from which quantity the error no longer becomes significantly smaller. Analogously, a sensitivity analysis should be performed regarding the sample position in the plate, i.e., which samples at which positions make the result particularly error-prone if they are omitted. Especially the idea of having less of the specimens at the borders of the plaque and more in a central area could be promising, particularly for the invariant decomposition method as the quantitative error seemed to be higher towards the edges of the plaque.
Another approach would be to deal with fourth-order tensors. In this work, only second level tensors were used, which already require a closure, which in principle means a reduction of information. The development or comparison of different methods for the interpolation of fourth level tensors would be highly interesting, but mathematically more complicated.
The ANN could be optimized in terms of the amount of layers, the depth of layers, the learning rate, the batch size, etc. However, as mentioned before, the authors strongly suggest an increase in training data to optimize the network for a specific process, in which similar fiber orientations are to be expected per plate/part instead of using it with this amount of data.
In the further course of the research project, the use of an interpolation method for three-dimensional components is envisaged. The ANN can be used for the 3D case without any changes, when inserting changing z values in the input file instead of constant ones. The component averaging method and the decomposition method can be used for 3D cases basically the same way as well. The only difference appears in the weighting. The inverse distance weighting method that was implemented for this paper, only reads in an and y value. The procedure can be similarly applied with just one more coordinate value. The rest of the methods can stay the same as the tensors are read in as 3 × 3 matrices with the possibility to change in any axis. For a more complex tool and component geometry, the extension of the interpolation functionality when dealing with different z-components of the positions of the FOT can then be tested. Alternatively, slicing the FOT through the thickness could provide for more input data and varying z coordinates. This approach has been chosen by Sabiston et al. [34] before.
Another interesting problem would be to investigate the (on average) increased isotropy of the tensors interpolated by the ANN. While this relationship seems logical for the component averaging method, it is so far not explainable for the neural network. A sensitivity analysis of both the ANN parameters and the input files and their effects on the interpolated tensors would be of interest in this context. Under certain circumstances, the additional, irrelevant input of a constant z-value for the 2D plate might be detrimental to the result.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The raw and processed data (μCT scans of all specimens and FOT of these specimens) as well as the Python scripts for all three methods required to reproduce the above findings have been published by the authors [44] and are available to download via the following DOI: https://doi.org/10.5445/IR/1000153725. This DOI includes all raw data acquired by the authors and all code written by the authors of this paper.
However, the up-to-date code for the component averaging and the decomposition method of second (and soon also fourth) order can also be found in the following github repository: https://github.com/ jewelsbla/oriopy. The methods have also been processed into a Python package named oriopy which can be accessed, downloaded and locally installed: https://pypi.org/project/oriopy/.
The C++ code to generate the FOT from the μCT scans is available to download from https://sourceforge.net/p/composight/code/HEAD/ tree/trunk/SiOTo/StructureTensorOrientation/FibreOrientation/Struct ureTensorOrientationFilter.cxx#l186. For further explanation please consider the paper by Pinter et al. cited in this work [35].
The Matlab method for the visualization of the tensors via tensor glyphs is available to download from https://de.mathworks.com/ma tlabcentral/fileexchange/27462-diffusion-tensor-field-dti-visualization. For more information please consider the paper by Barmpoutis et al. cited in this paper [39].