One of the main goals of the present paper is to describe a procedure for finding accurate parametric representations of real, experimental image data. To that end, a 3D microstructural mapping was carried out on a 1.4 mm-diameter cylinder of Al-5 wt%Cu, which was cut out of a cold-rolled plate (50% thickness reduction) that had been subsequently homogenized at 500°C for 24 h in air. The shape of individual grains in the specimen and the location of internal grain boundaries were revealed by 3DXRD microscopy measurements, performed at beamline BL20XU of the Japanese synchrotron radiation facility SPring-8 using a monochromatic beam of 32 keV x-rays (Poulsen, 2004). For 10 min prior to this room-temperature mapping, the specimen was subjected to a heat treatment at 575°C in air, which results in a liquid AlCu phase of approximately 2 vol% wetting the boundaries between the solid, aluminum-rich grains. Owing to the simultaneous presence of two phases, the resulting evolution of the sample’s microstructure is classified as Ostwald ripening (Wang and Glicksman, 2007). Once the sample is removed from the furnace, however, the liquid layer crystallizes and the growth/shrinkage of individual grains ceases.

Reconstruction of the 3DXRD data followed the protocol described in Dake et al. (2016), relying on the data processing routines of Schmidt (2005, 2014). To each voxel in the reconstructed volume, the software assigns the crystal lattice orientation that generates the most complete diffraction signal, whereby “completeness” is defined as the ratio between the number of experimentally detected diffraction spots associated with the voxel in question and the number of diffraction spots that are simulated to arise from this particular voxel if it were to have the assumed orientation. The grain labels were then assigned voxel-by-voxel to the orientation having the greatest completeness value. Formally, we describe the resulting image dataset as a mapping I E , r a w : W E , r a w = 1 , … , 531 × 1 , … , 321 × 1 , … , 321 → 0 , … , 943 , where each voxel coordinate is mapped to the corresponding grain label. Here, the label 0 is assigned to the background (i.e., voxels located outside the specimen). Each of the remaining labels is associated with one of the 943 grains.

However, with this reconstruction procedure the grain boundaries may manifest irregularities, such as local roughness, “island” voxels, zigzag shapes, or regions of fluctuating curvature (see Figure 1A) as a result of measurement uncertainties. These artifacts can be eliminated by treating the raw reconstruction as the initial configuration of a computational simulation of curvature-driven grain growth. If the duration of such a simulation is kept short enough, any boundary location manifesting severe curvature will tend to smoothen out, and any island voxels will be consumed by the surrounding grain, but no long-range translation of boundaries will occur—see Figure 1B. In the present paper, we employed 25 iterations of a 3D phase field algorithm (Krill and Chen, 2002) to reduce the roughness of grain boundaries in the raw 3DXRD reconstructions. The resulting smoothed experimental image is referred to as I E , s m o o t h : W E , s m o o t h → 0 , … , 938 , where W E , s m o o t h = W E , r a w . Note that some smaller grains vanished during the data smoothing procedure; consequently, I _E,smooth had fewer grains than I _E,raw (938 grains instead of 943).

In order to represent the 3D grain architecture of (measured and simulated) image data in an efficient way, we apply an optimization method to decompose the volume of interest into subvolumes using tessellations. Thus, to begin with, we briefly describe the tessellation models considered in the present paper. For additional details on tessellations in general, we refer, e.g., to Chiu et al. (2013).

Roughly speaking, a tessellation is a partitioning of space into pairwise disjoint sets, so-called cells. More precisely, a tessellation T in a sampling window W ⊂ R 3 is a countable collection of sets (cells), T = { C i T ⊂ W : i = 1,2 , … } , such that

1) int ( C i T ) ∩ int ( C j T ) = ∅  for  i ≠ j ,

For practical purposes, such as finding simplified representations of experimental image data, parametric tessellation models are probably the most suitable class of tessellations. The tessellation models considered in the present paper have in common that their cells are defined in terms of a distance function d T : R 3 × G → {R} , where { G } denotes the domain of generators (i.e., a set of admissible parameters of a single tessellation cell). For a (finite) set of generators G = { g i } i = 1 n T ⊂ G n T , the i-th cell C i T of a distance-based tessellation T is given by C i T = x ∈ W : d T x , g i ≤ d T x , g j  for each  j = 1 , … , n T . (1)

For brevity, we use the notation C i T x = 1 if  x ∈ C i T , 0 otherwise, to indicate whether a point x ∈ W belongs to the i-th cell, where i = 1 , … , n T .

The simplest model of a distance-based tessellation is the Voronoi tessellation, where d T ( x , s ) = ‖ x − s ‖ for x ∈ W with a generator s ∈ R 3 = G , and ‖x − s‖ denotes the Euclidean norm of x − s. While widely studied in literature, see e.g. Aurenhammer et al. (2013); Møller (1994); Okabe et al. (2000), the Voronoi tessellation is often found to be insufficiently flexible to fit experimental maps of polycrystalline materials (Šedivý et al., 2018); thus, more sophisticated tessellation models are needed. The fitting procedure considered in the present paper is able to handle many tessellations of the form given by Eq. 1 for which the distance function d T is differentiable. However, we focus on tessellation models that are special cases of generalized balanced power diagrams (GBPDs), listed here in order from simplest to most complex:

1) The Laguerre tessellation (Lautensack and Zuyev, 2008) is obtained if d T ( x , ( s , w ) ) = ‖ x − s ‖ 2 − w for x ∈ W , with a generator consisting of a seed point s ∈ R 3 and an additive weight w ∈ R .

Note that all of these models, except the Laguerre tessellation, can exhibit curved cell boundaries and thus non-convex cells. With this property comes the possibility, however, that cells are no longer connected, which might be undesirable when seeking parametric representations of polycrystalline materials. To rectify this issue, modifications of the original tessellation models, such as the one given by Šedivý et al. (2018), can be applied to fitted generators as a post-processing step. Another problem that affects all of the tessellation models described above is the possibility for a generator not to produce a corresponding cell. This can be mitigated by considering a volume-based cost function during the fitting that penalizes missing cells—see Section 2.3.

In this section we describe an efficient, gradient descent-based fitting procedure for GBPD-type tessellations. This procedure was originally introduced in Furat et al. (2021), but in Section 3 it will be applied to a broader class of tessellation models than in Furat et al. (2021).

Note that the fitting of a tessellation T = { C i T } i = 1 n T can be achieved by finding generators such that the similarity between the tessellation T and the ground truth image data is maximized. Formally, we consider the i-th grain of the ground truth image data as a map C i G T : Z 3 → { 0,1 } given by C i G T x = 1 if  x  belongs to the  i - th grain , 0 otherwise, with Z the set of all integers, i = 1, …, n _GT, and n _GT the number of grains in the sampling window W ⊂ R 3 . Nearest-neighbor interpolation can be used to extend the domain of C i G T from the integer lattice Z 3 to the continuous Euclidean space R 3 . Furthermore, let X F = { x j F } j = 1 n F ∈ Z 3 × n F be the set of all coordinates of voxels that belong to one of the grains, which we call the foreground voxels of the image data. If n _F denotes the number of foreground voxels in W , then for each j = 1, …, n _F there is an integer i = 1, …, n _GT such that C i G T ( x j F ) = 1 . In Section 3 we will consider the smoothed experimental image data from Section 2.1 (among others) and set C i G T x = 1 if  x ∈ W E , s m o o t h  and  I E , s m o o t h x = i , 0 otherwise, and X F = { x ∈ W E , s m o o t h : I E , s m o o t h ( x ) > 0 } .

Probably the most natural way to define the similarity between a tessellation T and the ground truth image data is to count the voxels at which each cell of the tessellation and the corresponding grain of the ground truth dataset overlap. To be more precise, the value of the objective function E : G n G T → [ 0 , ∞ ) for a set of generators G = { g i } i = 1 n G T with g i ∈ G is given by E G = 1 n F ∑ j = 1 n F ∑ i = 1 n G T C i T x j F C i G T x j F , (2) where the cells C 1 T , … , C n G T T of the tessellation T depend on the choice of the generators in G subject to n T = n G T . The corresponding fitting problem is thus to determine an optimal set of generators G o p t defined as G o p t = argmax G E G . (3)

It is easy to see that C j T ( x F ) with x F ∈ X F can be reformulated as C j T x F = argmin j ∗ d T x F , g 1 , … , d T x F , g n T , (4) where argmin j * is the j-th component of the n T -dimensional vector-valued argmin function, i.e., argmin ∗ : R n T → 0,1 n T , argmin j ∗ z = 1 if  z j ≤ z i  for all  i = 1 , … , n T , 0 otherwise, with z = ( z 1 , … , z n T ) ∈ R n T . In cases where the minimum is not unique, i.e., there are indices j 1 , j 2 ∈ { 1 , … , n T } with j ₁ ≠ j ₂ and z j 1 = z j 2 , only the component with the smallest index is set equal to 1. The function argmax^∗ is defined analogously.

Even though the distance function d T is differentiable (with respect to the generators), the fact that argmin^∗ in Eq. 4 does not have a derivative makes the objective function E defined in Eq. 2 non-differentiable. This leaves us having to resort to derivative-free optimization algorithms to solve Eq. 3, which in most cases converge slower than gradient descent methods (Audet and Hare, 2017). In order to increase efficiency, we slightly deviate from the original tessellation formulation by replacing the argmin^∗ function in Eq. 4 with a “softmin^∗” function—i.e., a softmax^∗ function with a negative argument, C ̃ j T ( x F ) = softmax j ∗ − d T ( x F , g 1 ) , … , − d T ( x F , g n T ) for x F ∈ X F . Here, the n T -dimensional function softmax ∗ : R n T → 0,1 n T , softmax j ∗ z = exp z j ∑ k = 1 n T ⁡ exp z k with z = ( z 1 , … , z n T ) ∈ R n T is a smooth version of the argmax^∗ function. So, instead of returning a vector the components of which are either 0 or 1, the softmax^∗ function is a vector-valued map, the components of which are continuous functions with values between 0 and 1. In fact, the output vector softmax^∗ z for some argument z ∈ R n T defines a discrete probability measure (i.e., the values of all components are between 0 and 1 and their sum is equal to 1), which assigns the highest probability to the index j if z _j ≥ z _i for all i ∈ { 1 , … , n T } . Here, the last property can be understood in the sense that the softmax^∗ function preserves the maximum of the input vector. The largest value of a component of the vector ( C ̃ 1 T ( x F ) , … , C ̃ n T T ( x F ) ) is therefore the one whose corresponding generator has the shortest tessellation distance to the given evaluation point x F .

The benefit of applying the softmax^∗ function instead of using the tessellation distances directly is found in the fact that the output vector of the softmax^∗ function is normalized, and thus for each evaluation point x F only the relative changes in the tessellation distances to the generators are considered, providing the same scale (i.e., values in [0, 1]) for all evaluation points. This trick is often used for multi-class classification problems in machine learning (Goodfellow et al., 2016). Furthermore, note that since softmax^∗ is a composition of differentiable functions and is itself therefore differentiable, the function C ̃ j T is differentiable, as well, for each j ∈ { 1 , … , n T } .

Because—in contrast to C j T , which is either 0 or 1— C ̃ j T assumes continuous values, it is necessary to adapt the objective function. Consequently, instead of E defined as in Eq. 2, we consider the function E ̃ : G n G T → ( − ∞ , 0 ] , where E ̃ G = 1 n F ∑ j = 1 n F ∑ i = 1 n G T λ N C E C ̃ i T x j F , C i G T x j F (5) with the negative (binary) cross-entropy loss function λ _NCE : [0,1]² → ( − ∞, 0] given by λ N C E ( y ̃ , y ) = y ⁡ log y ̃ + ( 1 − y ) log ( 1 − y ̃ ) . Note that the cross-entropy loss is often used in machine learning (Goodfellow et al., 2016) to compare the output of a classifier to ground truth data, which is basically the same purpose it serves here: If an evaluation point x F belongs to the i-th grain (i.e., C i G T ( x F ) = 1 ), C ̃ i T ( x F ) also needs to be close to 1 in order to maximize λ _NCE, and vice versa. Note that the modified objective function E ̃ is differentiable with respect to the generators. The fitted set of generators G ̃ o p t can then be obtained by solving G ̃ o p t = argmax G E ̃ G . (6)

In summary, we reformulated the original fitting problem, Eq. 3, into the differentiable version given in Eq. 6. This allows us to employ fast, gradient-based optimization algorithms, such as the one used in the present work: the stochastic gradient descent algorithm Adam (Kingma and Ba, 2015) (applied to the negative objective function). The optimization is stopped after a maximum of 25 iterations through (random permutations of) the dataset or if the objective function E ̃ defined in Eq. 5 does not increase by more than 10⁻⁴ in 3 iterations. Note that for the discretization of a fitted GBPD-type tessellation, we compute the (unique) cell labels using the classical definition in Eq. 4. The software implementation for the fitting and the discretization is based on Tensorflow (Abadi et al., 2015), which allows for highly parallel and even GPU-accelerated computations.

In addition to the procedure described in Section 2.3, we mention two additional methods from the literature for fitting tessellations to image data, to which we will refer below. We start with the procedure for Laguerre tessellations introduced in Spettl et al. (2016), which was used to acquire the initial parameter configuration in Section 3. Furthermore, in order to compare the results of our method described in Section 2.3, we also employed a different method for the fast fitting of GBPDs that was originally developed in Teferra and Rowenhorst (2018).

To evaluate the goodness of fit of the tessellation { C i T } i = 1 n T fitted to the foreground voxels X F = { x j F } j = 1 n F of the ground truth image data { C i G T } i = 1 n G T (with n G T = n T ), we consider various performance measures. The fraction of correctly assigned voxels is given by F c = 1 n F ∑ j = 1 n F ∑ i = 1 n G T C i T x j F C i G T x j F .

Similarly, the fraction of correctly assigned boundary voxels is defined as F c B = 1 n B ∑ j = 1 n B ∑ i = 1 n G T C i T x j B C i G T x j B , where X B = x j B j = 1 n B = { x 1 B ∈ X F : ‖ x 1 B − x 2 B ‖ ≤ 3  and  C i G T x 1 B = C ℓ G T x 2 B = 1 for some  x 2 B ∈ X F  and  i , ℓ = 1 , … , n G T , i ≠ ℓ } are the coordinates of the n _B grain boundary voxels (with respect to the 26-neighborhood, where the voxels x ₁, x ₂ are neighbors if ‖ x 1 − x 2 ‖ ≤ 3 ). Moreover, note that the fraction of empty cells F ₀ can be written as F 0 = 1 n T ∑ i = 1 n T 1 { C i T x F = 0  for all  x F ∈ X F } , where 1 { ⋅ } denotes the indicator function. Furthermore, consider the set of all grains that are a neighbor of the i-th grain in the ground truth image data (with respect to the 6-neighborhood of each voxel, where the voxels x ₁, x ₂ are neighbors if ‖x ₁ − x ₂‖ ≤ 1), N i G T = { ℓ ∈ 1 , … , n G T : C i G T x 1 F = C ℓ G T x 2 F = 1  for any  x 1 F , x 2 F ∈ X F  with  ‖ x 1 F − x 2 F ‖ ≤ 1 , i ≠ ℓ } , the set of all cells that are a neighbor of the i-th tessellation cell, N i T (defined analogously), and the resulting set of correctly assigned neighbors, N i = N i   G T ∩ N i T . Then, the fraction of cells for which all cell neighbors are correct can be written as N 0 = 1 n G T ∑ i = 1 n G T 1 { N i   G T = N i T } , and the mean number of incorrect cell neighbors is N ̄ = 1 n G T ∑ i = 1 n G T # N i   G T − # N i , where # denotes cardinality. These performance measures were calculated for simulated and experimental image data with both smooth and rough grain boundaries (see Sections 3.2, 3.3).

In this section, the fitting of tessellations to simulated image data is investigated. This allows us to study scenarios in which, in principle, the tessellations can perfectly describe the image data, which is usually not the case for experimental data. Apart from that, it is also possible to simulate the effect of noisy image data while still having access to the true grain boundaries.

In this section, the fitting method described in Section 2.3 is tested with realistic grain boundaries. For this purpose, experimental image data obtained from a 3DXRD mapping of an AlCu sample (Section 2.1) was used. As with the simulated data of Section 3.2, we first consider image data with smooth grain boundaries before tackling a dataset with rougher boundaries (Figure 1). The goal is to assess the robustness of the fitting procedure with respect to real-world grain boundary perturbations—originating, e.g., from measurement uncertainties—and also to determine whether the custom of preprocessing raw experimental image data to obtain smoother (and, thus, physically more sensible) grain boundaries prior to fitting is actually necessary.

As seen in Section 3.2.1, the fitted multiplicatively weighted Laguerre tessellation, the diagonal GBPD and the GBPD matched the smooth simulated image data very well, and a near perfect voxelwise accuracy was obtained. Between these tessellation models, there were only minor differences in the goodness of fit. Most notably, the GBPD was slightly worse at reconstructing the boundary voxels (see Table 1). Nevertheless, no significant decline in goodness of fit was observed even for the tessellation models that employ more parameters than necessary to describe the ground truth image data (which was generated from a multiplicatively weighted Laguerre tessellation). Furthermore, except for the number of empty cells, the gradient descent-based fitting procedure described in Section 2.3 achieved a notably better fit than the good results obtained by the direct approach of Teferra and Rowenhorst (2018) (see Section 2.4.2). On the other hand, in light of the results given in Table 1 for the initial Laguerre tessellation and the improved generators that resulted from the fitting procedure of Section 2.3, it is clear that the Laguerre tessellations with their flat boundaries (see Figure 3B) lacked sufficient flexibility for an accurate reconstruction of the ground truth data. However, the gradient descent-based fitting procedure still managed to bring about a significant improvement compared to the initial generators from the cross-entropy approach. This can likely be traced to the fact that the gradient descent-based approach considers a volume-based objective function (see Section 2.3) rather than an interface-based one (see Section 2.4.1). In general, we cannot expect to solve such a high-dimensional optimization problem by finding its global optimum; this would be equivalent to finding a perfect reconstruction of the multiplicatively weighted Laguerre tessellation that underlays the ground truth image data. Nevertheless, the fitted tessellations were quite close to the optimum, despite having been obtained by a local optimization method.

When it comes to the perturbed simulated image data investigated in Section 3.2.2, it is somewhat surprising how well the voxels of the smooth simulated image data could be reconstructed even from very noisy input image data (see Figure 7). As the same effect is observed for all three fitting methods considered in the present paper—i.e., the cross-entropy approach (Section 2.4.1), the direct approach (Section 2.4.2), and the gradient descent-based method (Section 2.3)—we attribute this robustness against perturbations to the inherent smoothing property of tessellations. Another observation that might surprise is the finding that the results for the direct approach were practically independent of the perturbation spread σ, but, just as in the case of the smooth simulated dataset, the gradient descent-based fitting procedure of Section 2.3 was able to surpass the direct method. The proposed method also achieved a better accuracy of the boundary voxels. In the latter case, however, some degradation could be observed with increasing σ. Some part of this degradation can be explained by the fact that a procedure producing more accurate approximations of the grain boundaries in the first place is going to be more sensitive to noise in the grain boundaries. Naturally, the deterioration in the fit of the boundary voxels also influences the considered grain neighborhood characteristics N ₀ and N ̄ —compare Tables 1, 2. Nevertheless, the decline in goodness of fit with σ is still well within reason, given the high noise level of the input image data (see Figure 6). In fact, from a visual comparison of the raw experimental image data in Figure 1B to the perturbed simulated image data in Figure 6, it can be seen that the case with the lowest level of noise (i.e., σ = 0) comes closest to the considered experimental data. This indicates that the gradient descent-based fitting procedure described in Section 2.3 is not only able to deal with the present level of measurement artifacts but also with much noisier scenarios.

For the experimental image data investigated in Sections 3.3.1, 3.3.2, quite high voxelwise accuracies F _c and F c B were observed, particularly for the more sophisticated tessellation models, as they are better at describing non-convex grain morphologies. Despite this fact, with respect to F _c and F c B the simpler Laguerre tessellation does a better job of fitting the experimental datasets than the simulated image—compare Tables 3, 4; Table 1. The same is true when considering the grain neighborhood characteristics N ₀ and N ̄ . As could already be anticipated from our analysis of the perturbed simulated image data, there was no significant difference in goodness of fit to the smoothed versus to the raw experimental datasets (for all performance measures). This leads us to conclude that there is no benefit to smoothing the grain boundaries in experimental image data prior to tessellation fitting.

One of the main advantages of the gradient descent-based fitting procedure described in Section 2.3 compared to other optimization approaches lies in its runtime performance. This comes from the reformulation of the objective function and the resulting ability to employ efficient gradient descent optimization algorithms. Another reason is the fact that the objective function and its gradient can both be computed on multiple CPU cores or even on GPUs. This opens up the possibility of reducing the (wall clock) time for the fitting procedure by employing hardware having a higher degree of parallelism (such as GPUs), which would not be so readily feasible for sequential fitting procedures. As runtime benchmarks are notorious for their dependence on a multitude of factors (in this case, for example, on the number of generators/grains, the number of voxels, the computer hardware, etc.), the runtimes quoted in Table 5 for the fitting of tessellations to the smoothed experimental image data should be taken with a grain of salt. In this particular case, the slower fitting of the GBPD model than the other tessellations can be attributed not only to the increased number of parameters but, more importantly, to the fact that a less efficient software implementation had to be used to compute the GBPD distance function. That being said, to achieve such a high goodness of fit for a dataset with 531 × 321 × 321 voxels and 938 grains, the runtimes are quite competitive, especially for the other tessellation models.

The fitting results for the direct approach of Teferra and Rowenhorst (2018) (Section 2.4.2) were discussed in the previous section. For all datasets, the direct approach delivered reasonably good fits of GBPDs, but they were consistently worse than those obtained by the gradient descent-based fitting procedure of Section 2.3. The only exception here is the fact that the direct method produced very few empty cells. Its main benefit, however, is the very short runtime of only a couple of seconds. Therefore, the direct approach is a good choice if a tessellation must be found quickly, but if the focus lies on the quality of fit, the gradient descent-based method developed in the present paper may be more suitable.

In Šedivý et al. (2016), still another method for fitting GBPDs was proposed, in which—similar to the present paper—a volume-based objective function is optimized. However, instead of a gradient descent method, Šedivý et al. employed a global stochastic optimization technique, the simulated annealing algorithm. As the name implies, this technique is inspired by the annealing (heat treatment) of metals. Specifically, during each iteration of the optimization, a random modification of the previous generators is proposed. These changes are accepted depending on whether the objective function is improved as well as on the current value of the ‘temperature.’ Here, the parameter ‘temperature’ governs the likelihood that a change in generators is accepted even though it leads to a worse value of the objective function. As the temperature is decreased over the course of the optimization, the probability increases that only improvements in the fit are accepted. In Šedivý et al. (2016), this fitting procedure was applied to both simulated and experimental image data. For the former case, which was a realization of a random GBPD, the quality of fit of the simulated annealing approach (F _c = 0.975, N ₀ = 0.604, N ̄ = 0.57 ) was similar to the results obtained in Section 3.2.1 of the present paper. However, since the simulated datasets were drawn from two different models (GBPD vs. multiplicatively weighted Laguerre tessellation), direct comparison warrants caution. Regarding the runtime of the fitting routine, it is mentioned in Šedivý et al. (2016) that roughly 19 h were needed to carry out 10 million optimization iterations on an Intel Xeon E3-1240 CPU with four 3.4 GHz cores (slightly slower than the Intel Core i7-4770K used in the present paper). The procedure stopped after 13 million iterations, which corresponds to a total runtime of about 24 h. When comparing this to the 16:42 h that the same fitting took in the present paper (see Table 5), we note that the dataset considered by Šedivý et al. had only about a 10th the size (180³ vs. 531 × 321 × 321 voxels) but more than twice as many grains (1894 vs. 938 grains). A direct comparison to their experimental dataset is omitted, as the two datasets are quite different. In summary, the goodness of fit achieved by the simulated annealing approach applied to simulated data was similar to that achieved by the gradient descent-based fitting procedure of the present paper, but the runtime reported in Šedivý et al. (2016) was significantly longer.