Optimized reconstruction of the crystallographic orientation density function based on a reduced set of orientations

In this work, a new method is developed to reconstruct the orientation distribution function from experimental data with a set of equally weighted orientations. It is based on the deterministic integer approximation method, but it minimizes the shortcomings of this method by optimizing the orientation grid and kernel function.

Optimized reconstruction of the crystallographic orientation density function based on a reduced set of orientations Abhishek Biswas, a * Napat Vajragupta, a Ralf Hielscher b and Alexander Hartmaier a a Interdisciplinary Centre for Advanced Materials Simulation, Ruhr University, 44801 Bochum, Germany, and b Chair of Applied Functional Analysis, Chemnitz University of Technology, 09126 Chemnitz, Germany. *Correspondence e-mail: abhishek.biswas@rub.de Crystallographic textures, as they develop for example during cold forming, can have a significant influence on the mechanical properties of metals, such as plastic anisotropy. Textures are typically characterized by a non-uniform distribution of crystallographic orientations that can be measured by diffraction experiments like electron backscatter diffraction (EBSD). Such experimental data usually contain a large number of data points, which must be significantly reduced to be used for numerical modeling. However, the challenge in such data reduction is to preserve the important characteristics of the experimental data, while reducing the volume and preserving the computational efficiency of the numerical model. For example, in micromechanical modeling, representative volume elements (RVEs) of the real microstructure are generated and the mechanical properties of these RVEs are studied by the crystal plasticity finite element method. In this work, a new method is developed for extracting a reduced set of orientations from EBSD data containing a large number of orientations. This approach is based on the established integer approximation method and it minimizes its shortcomings. Furthermore, the L 1 norm is applied as an error function; this is commonly used in texture analysis for quantitative assessment of the degree of approximation and can be used to control the convergence behavior. The method is tested on four experimental data sets to demonstrate its capabilities. This new method for the purposeful reduction of a set of orientations into equally weighted orientations is not only suitable for numerical simulation but also shows improvement in results in comparison with other available methods.

Introduction
Microstructure characterization of polycrystalline material performed using diffraction experiments provides information regarding the crystallographic orientations. Depending on the resolution and the size of the scanned area, a diffraction experiment like electron backscatter diffraction (EBSD) may consist of thousands and sometimes millions of measurement points of crystallographic orientations. A non-uniform distribution of these crystallographic orientations is represented with the help of a probability density function commonly referred to as the orientation density function (ODF) (Kocks et al., 2000). The ODF can then be estimated by superposition of radially symmetric kernel functions centered at each crystallographic orientation (Hielscher, 2013).
Texture plays an important role in material behavior and is mainly responsible for anisotropy (Hielscher et al., 2014). Numerical predictions of the anisotropic mechanical response in metals have been of special interest. In view of this, many macroscopic numerical models have been proposed [a few examples of such models are given by Hill (1948) and Barlat et al. (1991)]. However, these models lack a clear connection to the important microstructural features like texture, grain size etc. that influence the mechanical response. Micromechanical modeling approaches address these limitations by incorporating the important microstructural features into the numerical model.
To properly characterize the crystallographic orientations of grains in a polycrystalline material, EBSD measurements are usually performed with high spatial resolution (Humphreys, 2004). This inevitably results in a large number of measurement points. For micromechanical modeling of polycrystalline materials in the crystal plasticity finite element method (CPFEM) framework (Roters et al., 2010), an optimal number of orientations that closely represent the material characteristics are required. There are two main approaches for modeling polycrystals. In the first approach the grain volume fractions are considered equal and fixed, which requires a set of equally weighted orientations. In this regard, various ODF reconstruction strategies have been proposed. One such method is the 'hybrid integer approximation', which has been introduced for discrete ODFs by Eisenlohr & Roters (2008). This method combines the probabilistic sampling approach of Tó th & Houtte (1992) with the deterministic integer approximation of Leffers & Jensen (1986).
In the second modeling approach, grain volume statistical distributions are considered, and thus the reconstruction of the ODF is done by a set of weighted orientations. Recent work by Schaeben et al. (2017) introduces such a method, which utilizes the Dirichlet kernel to provide an unbiased estimate of the Fourier coefficient up to any finite order. The corresponding weights are estimated numerically and can be related to the Fourier transform due to the linearity of both the spatial and spectral domains. Further information about other reconstruction methods can be found in the work of Eisenlohr & Roters (2008), Schaeben et al. (2017) and references therein.
In the following, we introduce a method for reconstructing the ODF estimated from a given set of experimentally measured orientations by a smaller number of equally weighted orientations (sample orientations). The method begins with an initial guess of an equi-spaced SO(3) grid. The experimentally measured orientations are superimposed on this grid. The orientations lying within a distance of half grid spacing from a given grid point are defined as a cluster around that grid point. By utilizing an iterative scheme based on the work of Leffers & Jensen (1986), the number of required sample orientations is divided among the grid points in proportion to the size of the cluster of input orientations and it is represented by an integer value. Any cluster of orientations around a grid point with an associated integer value greater than one is further divided into smaller groups of orientations from which mean orientations are sampled. These mean orientations collectively form the reduced set of orientations. The ODF is estimated from the reduced set of orientations by employing a kernel function with an optimized shape parameter. The error function is the L 1 norm of the difference between the ODFs from the experiments and the reduced set of orientations. This process is continued iteratively with the intention of minimizing the error function by optimizing the SO(3) grid. This algorithm is implemented using the MTEX (Bachmann et al., 2010) toolbox in MATLAB.
The capability of this method is demonstrated by implementing it on three cubic crystal cases in the form of coldrolled OFHC copper (Rolled-Cu) (Anand, 2004), additively manufactured 316L stainless steel produced by two different laser powers -1000 W (316L-1000W) and 400 W (316L-400W) (Biswas et al., 2019) -and one orthorhombic case, forsterite, which is available as an EBSD example in MTEX (Bachmann et al., 2010). The orientations for Rolled-Cu are generated by plane strain compression simulation as suggested in the work of Anand (2004) and by extracting the grain orientations at every integration point in the last time step. The simulation is then repeated five times with different sets of random input orientations and the results from all simulations are combined to obtain 27 000 orientations.

Problem setup
Estimation of an ODF f : SOð3Þ ! R from crystallographic orientation measurements g i , i ¼ 1; . . . ; N, is a classical problem in crystallographic texture analysis (see Table 1 for notation). The most common method is called kernel density estimation (cf. Hielscher, 2013) and uses a kernel function : ½0; ! R to estimate an ODF as where ffðg i ; gÞ denotes the disorientation angle between the orientations g i and g. Note that the estimated ODF f heavily depends on the choice of the kernel function . Here we Table 1 Inventory of symbols and notations. Number of reduced orientations that fall into each of the cells q k , k ¼ 1; . . . ; K Equally spaced grid of orientations used for estimating the error between two ODFs restrict ourselves to the de la Vallé e Poussin kernel (cf. Schaeben, 1997), the half-width of which can be controlled by the parameter > 1. The objective in this work is to find for a given N a much smaller numberÑ N ( N of orientationsg g i ; i ¼ 1; . . . ;Ñ N, and a kernel parameter such that the ODF is a good approximation of the ODF f estimated from the full data set, i.e.f f ' f . As an error measure we consider the L 1 norm, which has been used previously in texture analysis by Schaeben et al. (2017) and Bozzolo et al. (2007).
measures the volume of differently oriented orientations between the initially estimated ODF f and its approximationf f . In the extreme case when the ODFs f andf f are concentrated around disjoint orientations the error approaches 2. On the other hand, whenf f approximates f better the error approaches zero.

Outline of the algorithm
The idea of our algorithm can be summarized in the following steps: (1) Estimate an initial ODF f from the input orientations g i , i ¼ 1; . . . ; N, using a fixed kernel function (Section 2.1).
(2) Fix a subdivision of the orientation space into M cells with resolution Á (Section 2.3).
(4) Find a scaling factor (0 < < 1) such that the integersñ n j ¼ roundð Â n j Þ satisfy the condition P M j¼1ñ n j 'Ñ N, whereÑ N is the target number of orientations of the reduced set (Section 2.4).
(7) Optimize the kernel function shape parameter such that the misfit k f Àf f k 1 is minimized (Section 2.6).
(8) Repeat steps 2 to 7 for different subdivisions of the orientation space and find the optimal resolution Á with respect to the final misfit k f Àf f k 1 (Section 2.7).

Subdivision of the input orientations
The process begins with fixing a resolution Á and a corresponding decomposition of the Euler angle space into M cells of approximately equal volumes corresponding to each grid point (cf. Fig. 1). Next, we count the number n j , j ¼ 1; . . . ; M, of orientations that fall into each of these cells.
An illustration of this process is shown in Fig. 1 for the test case of Rolled-Cu orientations. For ease of visualization, a grid spacing of Á = 15 is selected. The equi-spaced SO (3) is represented by black dots and the input orientations by red dots.

Integer approximation
In the second step of this algorithm, the target number of samplesÑ N is distributed among the cells in proportion to their counts n j , j ¼ 1; . . . ; M. This step is based on the work of Eisenlohr & Roters (2008), which utilized the original idea of Leffers & Jensen (1986). It is referred to as 'integer approximation (IA)' in the work of Eisenlohr & Roters (2008). Henceforth, we will use this term accordingly. This iterative scheme estimates a scaling factor , using a binary search algorithm as suggested by Eisenlohr & Roters (2008), which is then multiplied with the counts n j and rounded to the closest integer:ñ subject to the condition that

Determination of the downsampling
The IA provides an integer numberñ n j of orientations that need to be computed for each subdivision cell. The objective of this step is to reduce the n j orientations in each cell to onlỹ n n j orientations such that the mean of the new orientations is similar to the mean of the original orientations. Hence, if n n j ¼ 1, we simply take the mean of all orientations within the cell as the new orientationg g j .
For cases withñ n j > 1, we subdivide the original orientations g j i , i ¼ 1; . . . ; n j , within each cell j, j ¼ 1; . . . ; M, intoñ n j clusters and compute the new orientationsg g j 1 ; . . . ;g g jñ n j as the mean orientations of each of the clusters. The subdivision is performed randomly; other non-random subdivision methods like sequential sampling were also tested for all of the test cases, but the error k f Àf f k 1 in all cases varied within 1%. To keep the subdivision process unbiased, the random method is chosen. This is visualized in Fig. 1, in which a grid point is shown with a cluster of orientations (red dots). In this case, there are 117 orientations in the cell; if we suppose that the correspondingñ n j ¼ 4, then these orientations are randomly subdivided into four approximately equal clusters containing 29, 30, 29, 29 orientations. The 117 original orientations are then approximated by theñ n j ¼ 4 mean orientations of each cluster.
Let us further illustrate this process on a 1D example. Here the input data consist of 200 points ranging from 0 to 1. The corresponding density function f is visualized in Fig. 2 (plotted in black). The input data are subdivided into 20 equally spaced bins and the number of data points corresponding to each bin is fed to the IA method which reduces the data set to only 25 points. The density function estimated from this downsampled data set is plotted in green in Fig. 2. At this stage k f Àf f k 1 ¼ 0:29 for the 1D example. This 1D example is used in the later sections for the purpose of illustration.
2.6. Optimal shape parameter of the kernel function In the above 1D example we observe that the density function recovered from the reduced data set is much too oscillatory. This is because we have chosen the kernel functioñ shape parameter for the reduced data set to be equal to the shape parameter of the kernel function that we used for the original data set. However, in general, a smaller sample size requires a broader kernel function. Hence, the objective of this section is to optimize such that the corresponding ODF,f fits best the original ODF, The optimization process of starts with the initial guess equal to and progresses in fixed steps selecting the least value of k f Àf f k 1 . The effect of this optimization for the 1D example is shown in Fig. 3, in which k f Àf f k 1 reduces from 0.29 (Fig. 2)  Kernel density estimate (KDE) of the 1D example and the reconstructed data with 25 points using a Gaussian kernel. Results from various stages of optimization are shown as stems on the x axis which are sampled data (Ñ N ¼ 25) from various stages of optimization. Initially from only integer approximation k f Àf f k 1 ¼ 0:29 shown in green, after kernel optimization k f Àf f k 1 ¼ 0:17 in blue and finally after subdivision optimization k f Àf f k 1 ¼ 0:12 in red.

Figure 3
Effect of kernel shape parameter on k f Àf f k 1 between the input ( f ) and output (f f ) kernel density estimates for the 1D example. This comparison justifies the kernel shape optimization in the proposed algorithm.
The results of kernel optimization for various test cases and differentÑ N are shown in Table 2. As explained earlier, these results show that asÑ N increases the kernel function becomes finer (indicated by smaller values of ).

Optimum grid spacing
In the previous sections, we have found a subsampling of our input orientations for a given subdivision of the orientation space into cells and an optimal kernel function to recover an ODFf f. In this section, we aim to optimize the subdivision of the orientation space to obtain an optimal fit between the ODF f computed from the input orientations and the ODFf f computed from the reduced set of orientations. Fig. 4 indicates the dependency of k f Àf f k 1 on Á for the 1D problem, and the optimization of Á to minimize k f Àf f k 1 further yields a result of 0.12 shown in Fig. 2 (plotted as a red line). The optimization procedure is similar to the kernel optimization performed in Section 2.6. The algorithm proceeds in steps of 1 unit of Á and selects the least value of k f Àf f k 1 estimated over a range of Á.
Similarly, the effect of Á on k f Àf f k 1 for variousÑ N is shown in Fig. 5 for an actual EBSD data set (test case 316L-1000W).
We observe that the error increases in both directions, i.e. if the subdivision becomes too fine or if the subdivision becomes too coarse. On the one hand, if the subdivision becomes too coarse orientations within one cell may have a large misorientation angle and replacing them by their mean value increases the error. On the other hand, if the subdivision is too fine, the integer approximation will result in severe rounding errors.

Results
In this section the results of the reconstruction process (Table 3) are shown in the form of contour plots (Figs. 6-9) for all of the test cases. The results in Table 3 show that k f Àf f k 1 reduces with risingÑ N value; however, the corresponding computational time also increases. In these figures, the contour plots of the ODF from input orientations and the reconstructed ODF (forÑ N = 150, 400, 1000) are shown. In these contour plots, the peaks and contours in the input ODF are successfully captured by the reconstructed ODF, closely maintaining the intensity of the contours.

Figure 4
Effect of grid spacing Á on k f Àf f k 1 between the input ( f ) and output (f f ) kernel density estimates for the 1D example.   For a more detailed comparison between the ODFs f andf f the power plots estimated from the ODFs of the test cases are shown in Fig. 11. These power plots are estimated by summing the squared Fourier coefficients of a given harmonic order, which shows the contribution of each harmonic order to the texture index; a detailed mathematical description is given in the work of Schaeben et al. (2017). The harmonic contribution from the ODFf f closely matches the value from f. This is also observed in Fig. 10 in which the error k f Àf f k 1 reduces asÑ N increases.

Application
One of the main applications of this method is numerical modeling like micromechanical modeling. Here, an example of this process is presented for Rolled-Cu. The input is in the form of 27 000 crystallographic orientations. We choose a local crystal plasticity (CP) model without the effect of the strain gradient as described by Ma & Hartmaier (2014) for numerical modeling of material behavior. The material is assumed to be constructed of periodically repeating volume elements known as the representative volume element (RVE) (refer to Fig. 12). For details of the applied periodic boundary conditions and homogenization scheme, please refer to Vajragupta et al. (2017).
A virtual uniaxial test is performed for all the RVEs by applying displacement along the z or 33 direction. The CP parameters are fitted by comparing the homogenized virtual uniaxial tensile test from the RVE consisting of 27 000 grains with experimental data; therefore the entire input orientation set was used for this test, and no reconstruction process is used for this RVE. Keeping the CP parameters the same, smaller RVEs comprising 64, 216, 512, 1000, 3375, 4913, 8000 and 10 648 grains are generated and corresponding grain crystallographic orientations are generated with the reconstruction algorithm (Table 4) and hybrid IA method. The results from the smaller RVEs created using both the methods are compared with the reference RVE, i.e. the RVE with 27 000 grains.
To exclude the influence of the grain boundary misorientation on the mechanical response of the material, the Comparison of k f Àf f k 1 between input ( f ) and output (f f ) ODFs for the test cases: Rolled-Cu, 316L-1000W, 316L-400W and forsterite calculated for a number of extracted samples.   Biswas et al. (2019) was implemented by using the extracted samples and the RVE geometry; the target distribution followed can be found in the work of Mackenzie (1964).
Since these orientation sets reconstruct the same input ODF, the output from the virtual tensile test data should be comparable to the results obtained with the model consisting of 27 000 orientations. The finite element method (FEM) simulations are performed with ABAQUS (Simulia, 2012). Each grain in the RVE is discretized using a single C3D8 hexahedral element (Fig. 12).
The homogenized true stress-strain plot shown in Fig. 13 indicates the difference between the results of the RVEs consisting of discrete orientations generated by the reconstruction algorithm and the hybrid IA method. This mismatch can be attributed to the different input data set. In the case of the hybrid IA, the input is the discrete ODF calculated on predefined SO (3)  Representative volume elements withÑ N ¼ [64,216,512,1000,3375,8000] grains. Each grain is represented by one cube corresponding to one finite element.

Figure 13
Homogenized CPFEM result comparison for the case of Rolled-Cu for a few selected RVEs constructed using orientations from the reconstruction algorithm and the hybrid IA method; experimental data are obtained by digitizing Fig. 13(a) of Knezevic & Landry (2015). the deterministic IA, whereas the proposed reconstruction algorithm starts with the experimental data set.
Referring to the generated RVEs, the results indicate that the RVE with 216 grains is sufficient to predict the homogenized mechanical properties; a detailed examination of the local stress values (calculated at integration points) gives a different outlook in comparison with the homogenized mechanical properties. The local stress values can be important for modeling phenomena like fracture, damage etc.
The distribution of the local stress values in smaller RVEs is shown in Fig. 14 and compared with the local stress distribution in the reference RVE, and the difference is illustrated in Fig. 15 as the root-mean-square difference. Evidently, as the number of grains in the RVE increases, the difference between the stress distribution from the smaller RVE and the reference RVE also decreases. These results indicate that the value k f Àf f k 1 only specifies a minimum requirement, and the strategy for numerical modeling should be decided by the objective of modeling. If the objective in this test case (Rolled-Cu) is to model damage, perhaps an RVE with a minimum of 1000 grains should be considered.
Another important aspect of textured materials is the anisotropy in material properties. This study is performed for a hypothetical test case of copper with fiber texture in which the 110 fiber is parallel to the Y direction of the RVE. The CPFEM parameters are the same as those in the work of Anand (2004)  Comparison of stress components along the tensile load direction ( 33 ) (distribution in the form of a histogram) for RVEs with 64, 216, 512, 1000, 3375, 4913, 8000 and 10 648 grains with the reference RVE (27 000 grains). The bin width of the histogram is 10 MPa.

Figure 15
Root-mean-square difference between the stress distribution from the reference RVE (27 000  are similar to those shown in Fig. 12. The boundary conditions are similar to the previous test case of Rolled-Cu. However, in this test case, the uniaxial tensile test is performed in all three sample directions for all the RVEs to study the anisotropy in homogenized stress-strain behavior. The input orientation consists of 500 000 orientations generated analytically using MTEX (Bachmann et al., 2010); the results of the reconstruction process are shown in Table 5. Fig. 16 shows the comparison between the stress-strain plot of the RVEs. Since the 110 fiber is aligned along the Y direction the stress values are much higher than in the other two directions. All the RVEs show a similar stress-strain behavior. This indicates that the extracted samples can successfully reconstruct the input ODF and are also able to capture the mechanical behavior.

Conclusion
In this work, an algorithm for the reconstruction of the ODF from an EBSD experiment by a set of equally weighted orientations has been proposed. It is based on the deterministic integer approximation method introduced by Leffers & Jensen (1986), but the previously reported problem of overweighting is tackled by optimizing the SO(3) grid and kernel function used for the reconstruction. The quality of the reconstruction is judged not only by the L 1 norm of the difference between the input and the reconstructed ODF but also by the ODF power plot. The application of this method in the prediction of mechanical behavior using CPFEM provides further insight into the importance of precise representation of the reduced orientation set in CPFEM simulations. This study shows that the L 1 norm of the difference between the input and the reconstructed ODF estimated from the sample orientations gives only a minimum criterion for the number of samples to be incorporated in the RVE. However, the influence of sample size on the local mechanical output like stress can be observed and should be considered during micromechanical modeling. In addition to this, we also demonstrate the ability of the reduced orientation set to predict anisotropy in yield strength and hardening behavior, through micromechanical simulations of the RVE with fiber texture. This case study shows a good agreement for the prediction of these mechanical behaviors between the various sample sizes ranging from 216 to 10 648 orientations.

Figure 16
Homogenized CPFEM result comparison for the case of copper with fiber texture.