{'fraction_non_alphanumeric': 0.03444741057127603, 'fraction_numerical': 0.01025093432995195, 'mean_word_length': 4.664207088423854, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fitting Laguerre tessellation approximations to tomographic image data

The analysis of polycrystalline materials benefits greatly from accurate quantitative descriptions of their grain structures. Laguerre tessellations approximate such grain structures very well. However, it is a quite challenging problem to fit a Laguerre tessellation to tomographic data, as a high-dimensional optimization problem with many local minima must be solved. In this paper, we formulate a version of this optimization problem that can be solved quickly using the cross-entropy method, a robust stochastic optimization technique that can avoid becoming trapped in local minima. We demonstrate the effectiveness of our approach by applying it to both artificially generated and experimentally produced tomographic data.


Introduction
In recent years there have been significant advances in the tomographic characterization of materials. As a result, it is now possible to carry out detailed investigations of the 3D grain structures of polycrystalline materials; see, e.g., [1][2][3][4][5][6][7]. A fundamental ingredient in any such investigation is a suitable quantitative description of the grain morphology. Such a description contains the key features of the structure, ideally free from noise and imaging artifacts. A good description usually results in significant data compression, describing large 3D voxel data sets using only a small number of parameters. Data compression is necessary, for example, when carrying out analysis of sequences of tomographic data sets (e.g., the high time resolution in-situ synchrotron images considered in [8]). In addition, the description of tomographic data via tessellations provides a basis for the statistical analysis of grain structures and, in some cases, can be used to develop stochastic models of material microstructures; see, e.g., [8][9][10][11][12][13][14][15][16].
The most commonly used quantitative descriptions of space-filling grain ensembles are based on tessellations, which divide the space into disjoint regions called cells.
The cells represent the individual grains. The most widely used tessellation model is the Voronoi tessellation (see, e.g., [17,18]), which takes, as parameters, a set of generating points. The space is then divided into convex cells by assigning each point to its nearest generator. The Laguerre tessellation (see, e.g., [17,18]) is a generalization of the Voronoi tessellation that also partitions the space into convex cells. Like the Voronoi tessellation, the Laguerre tessellation is generated by a set of points; however, unlike the Voronoi tessellation, these points are weighted, with the weights influencing the size of the cells. Consequently, the Laguerre tessellation is able to describe a wider range of structures than the Voronoi tessellation. For this reason, the Laguerre tessellation is a popular choice for modeling polycrystalline grain structures [8,9,[19][20][21][22] and other materials, such as foams [12,23,24].
In order to describe a tessellation by a set of generating points, it is necessary to solve an inverse problem: that is, a set of generating points that produce the observed cells must be found. The Voronoi inverse problem (VIP) is well-studied; see, for example, [25][26][27][28][29][30]. Recently, Duan et al. [31] proposed an algorithm that finds solutions to the Laguerre inverse problem (LIP). Although the examples considered in [31] are restricted to 2D, the methodology is easily applied in higher dimensions.
The solutions to the VIP and the LIP assume that the empirical data constitute perfect descriptions of the observed cells. However, this is not true when working with tomographic data, which is distorted by noise and also contains imprecision arising from discretization during the imaging process. It is also worth noting that real-world materials are not perfectly described by Laguerre tessellations (though the descriptions can be quite good). These limitations mean that methods that attempt to invert a tessellation extracted from the tomographic data do not, in general, result in good fits. The LIP is solved by iteratively finding the generating points of the given tessellations. When applied to imperfect data, this iterative procedure propagates errors, resulting in tessellations that have little correspondence to the tomographic data. Thus, when dealing with empirical data, it is not appropriate to attempt to solve the LIP. Instead, the generating points of a Laguerre tessellation that is a good approximation of the material must be found. This is, at its core, an optimization problem. We call this problem the Laguerre approximation problem (LAP). The corresponding Voronoi approximation problem has been considered in the literature, beginning with [32].
A simple heuristic approach for solving the LAP was proposed in [6]. More sophisticated approaches, which formulate and solve an optimization problem, are described in [33][34][35]. Although these techniques provide good fits in certain settings, they are either limited to small sample sizes or require the considered tessellations to be sufficiently regular.
In this paper, we present a fast and robust algorithm for fitting Laguerre approximations to large data sets. More precisely, we formulate an optimization problem where we minimize the discrepancy between the grain boundaries observed in the image data and the grain boundaries produced by our Laguerre approximation. The cost function is chosen so that it can be evaluated very efficiently and that all necessary information can be easily obtained from image data. We then solve the optimization problem using the cross-entropy (CE) method [36][37][38][39], a stochastic optimization algorithm that is able to escape local minima. We carry out experiments on both real and artificially-generated image data that show our approach is able to produce very good fits.
This paper is structured as follows. In Section 2, we review some key properties of Laguerre tessellations. In Section 3, we give a more complete description of the LAP and formulate our optimization problem. Then, in Section 4, we introduce the CE method as a robust tool for solving this optimization problem. Section 5 gives results for both artificial and experimental data that demonstrate the effectiveness of our approach. Finally, Section 6 summarizes our results and suggests directions for further research.

The Laguerre tessellation
In the following section, we define Voronoi and Laguerre tessellations and give a number of properties that we will use to solve the LAP. For notational convenience, we only consider tessellations in R 3 . However, our methods are easily applied in other dimensions.
The Voronoi tessellation is defined by a locally finite set of generating points, {x i } i∈I , where I ⊆ N denotes the index set consisting of natural numbers. The cell corresponding to the ith generating point, x i ∈ R 3 , is given by where · is the Euclidean norm on R 3 . In the case of the Laguerre tessellation, the generating points, {(x i , r i )} i∈I , are weighted, where we assume that r i > 0. The cells of the Laguerre tessellation, {C i } i∈I , are then defined by for all i ∈ I, where the Euclidean norm used in (1) is replaced by the so-called power distance pow(y, (x, r)) = y − x 2 − r 2 .
Given the generator points, the faces of the cells can be computed efficiently; see, e.g., [40][41][42]. As we take the weights r i to be positive real numbers, the generating points have a geometric interpretation as spheres. That is, we can represent a generating point, (x i , r i ), as a sphere with center x i and radius r i . The flexibility of the Laguerre tessellation comes at a cost. For example, while each cell of a Voronoi tessellation contains its generating point, the generating points of a Laguerre tessellation may not be contained in their corresponding cells. In some cases, a generating point of a Laguerre tessellation may not even produce a cell. That is, it is possible that there exists a point, (x i , r i ), such that C i = ∅; see, for example, [40]. In addition, while the cells of a Voronoi tessellation uniquely determine its generating points, there are uncountably many sets of generating points that can generate a given Laguerre tessellation; see, [31,43,44]. Thus, while the VIP has a unique solution, the LIP has uncountably many solutions. Under mild conditions, however, it can be shown that the generating points of a Laguerre tessellation are uniquely determined given one generating point and the weight (or one coordinate) of a generating point in an adjacent cell; see [31].
We will often find it convenient to consider planes that are equidistant from two generating points (under their respective power distances). That is, given two generating points, (x i , r i ) and (x j , r j ), we consider the separating plane given by is the unit normal vector of P i,j and is the shortest distance from the origin to the plane; see, e.g., [43, Section 2.1] for more details. The plane P i,j defines a half-space which covers the cell C i . Note that the intersection of all such half-spaces, k∈I\{i} H i,k , defines the cell C i . Because the normal vectors are taken to be unit vectors, the distance from an arbitrary point, x ∈ R 3 , to the plane P i,j is given by

The Laguerre approximation problem (LAP)
Given an exact description of a Laguerre tessellation (e.g., in terms of the half-spaces defined above), it is not too difficult to solve the LIP. That is, it is straightforward to find a set of weighted points that is able to generate the given tessellation. Furthermore, these points can be chosen to satisfy certain constraints; see [31]. When dealing with tomographic data, however, the description of the tessellation is not exact. This is because, even if the material itself can be perfectly described by a Laguerre tessellation, the noise and discretization errors inherent in the imaging process mean that the cell boundaries extracted from the data are subject to error. The LIP can be solved as described in [31]. However, the iterative computation of generators is sensitive to imperfect data and errors propagate quickly. As a result, the ensuing tessellations often do not correspond well to the tomographic data. Therefore, solving the LIP for empirical data is generally ill advised. Instead, we wish to find generating points that produce a Laguerre tessellation that is as close as possible to the tomographic data (with respect to some metric). This is, at its core, an optimization problem. In order to properly formulate this optimization problem, a discrepancy measure needs to be defined. We then choose the generating points of the approximating tessellation in order to minimize this discrepancy. The choice of discrepancy measure depends on the nature of empirical data. In this paper, we work with tomographic data.

Tomographic image data
We assume that the tomographic image data constitutes a collection of voxels in a convex window. Furthermore, we assume that the image has already been segmented -e.g., by the watershed transform; see [45,46] -and that the voxels have been labeled by their corresponding grains. Thus, the data is of the form {I(x, y, z) ∈ {0, . . . , N } : (x, y, z) ∈ W }, where N ≥ 1 denotes the number of grains and W ⊂ N 3 is a grid of voxel coordinates. The ith grain region, R i , is then given by R i = {(x, y, z) ∈ W : I(x, y, z) = i}. Note that R 0 does not correspond to an actual grain, but is either the empty set or a collection of one voxel thick layers that separate grains. Such thin layers often arise in segmentation procedures such as the watershed transformation. A grain region R i may correspond to the grain itself or a region that contains the grain (if the space is not completely filled with grains), cf. Figure 1. We assume that the grain regions are roughly convex and that the segmentation is of a high quality.

Discrepancy measures
The most direct way to measure the discrepancy between 3D image data and a Laguerre tessellation generated by the points {(x i , r i )} i∈I is by counting the number Figure 1. Illustration of voxels labeled by their corresponding grains. The grain regions are assumed to form a tessellation with convex cells. The grains themselves may not fill their corresponding regions completely, as depicted for grain number 3. Here, we show the grain region (i.e., voxels with label 3), a cell approximating the grain region (colored cyan) and the grain itself (colored blue).
of incorrectly assigned voxels. That is, we calculate where I is a discretized version of the tessellation generated by {(x i , r i )} i∈I and #A denotes the number of elements in some set A. Many other discrepancy measures considered in the literature are of a similar form. For example, in [35], the difference between a cell in the approximating tessellation and its equivalent in the empirical data is measured by the intersection of the approximating cell with the corresponding adjacent cells in the empirical data. Minimizing the total overlap of the approximating cells is equivalent to minimizing the number of incorrectly assigned voxels. We call discrepancy measures that aim to minimize the volume of the difference between the empirical and approximating tessellations volume-based measures.
Although minimizing a volume-based measure gives a very good fit to the data, such measures are poorly suited to most numerical optimization methods. This is because evaluating such a discrepancy measure is computationally expensive. For example, when considering the number of incorrectly assigned voxels, a discretized version of the approximating tessellation needs to be generated. In practice, this limits the number of times such a discrepancy measure can be evaluated. As a result, approaches that aim to minimize volume-based discrepancies are restricted to small data sets with small numbers of grains (e.g., 109 grains in [34]) or are forced to use non-optimal optimization techniques, such as gradient-descent, which do not require too many evaluations of the discrepancy and, as such, limit the time taken by the fitting procedure (ideally to substantially less than 24 hours).
As we propose to use stochastic optimization techniques to solve this problem, we need a discrepancy measure that is significantly faster to evaluate. We can achieve this by considering an interface-based discrepancy measure -a measure that considers only the boundaries between cells -instead of a volume-based one. If an approximating tessellation can accurately reproduce these interfaces, it will be a good approximation of the data. The primary advantage of interface-based discrepancy measures is that they can be calculated from the generating points of the approximating tessellation without the need to generate the tessellation itself. In order to define such a discrepancy measure, we consider the sets of voxels that separate adjacent cells. We define the interface between two grains, i and j, in the empirical data by where N 26 (x, y, z) denotes the 26-neighborhood of (x, y, z) -i.e., the voxels (x , y , z ) that have a distance less than or equal to √ 3 from (x, y, z). This set contains all voxels that touch both grains. Note that N i,j = ∅ if the grains are not adjacent. If the tessellation generated by {(x i , r i )} i∈I is a good approximation to the empirical data, then the plane separating the generating points of adjacent cells i and j (which defines the boundary of the two cells) should be close to N i,j . Thus, we measure the distance between N i,j and P i,j , the plane separating cells i and j in the approximating tessellation. The discrepancy between the approximating tessellation and the empirical data is then given by the sum of squares of these distances. That is, we define where d(y, P j,k ) is given in (2) and {P j,k } j,k∈I is the set of separating planes determined by the generating points {(x i , r i )} i∈I . Note that this discrepancy measure can be calculated without generating the approximating tessellation. A separating plane can be computed for all pairs of generators (or generator candidates, when performing the optimization) -even if their cells are empty or not adjacent in the Laguerre tessellation. Therefore, such 'degenerate' cases are not a problem for our approach. By matching all separating planes to the corresponding test points at once, Laguerre cells computed based on 'good' configurations of generators will match the desired cell configuration well enough.

Minimizing the discrepancy
When aiming to minimize a volume-based discrepancy, the LAP reduces to the optimization problem of finding There are a number of significant difficulties which must be overcome in solving (4). In particular, (i) the optimization problem is high-dimensional, (ii) the optimization problem has many local minima, (iii) the discrepancy is expensive to evaluate.
The approach developed in [34] avoids problems (ii) and (iii) by finding the exact solution of a linear program (an optimization problem where a linear cost function is minimized subject to linear constraints). However, when the LAP is transformed into a linear program, the size of the resulting problem grows very quickly in both the number of voxels and the number of grains. This means that the linear programming approach cannot be applied to large 3D data sets containing many grains.
In [35], gradient-descent methods are used to obtain Laguerre approximations (with a relatively low number of evaluations of the discrepancy). Unlike the linear programming approach, the dimensionality of the optimization problem grows only linearly in the number of grains (and does not depend on the number of voxels). Thus, [35] avoids problem (iii) and, to a lesser extent, problem (i). However, the gradient-descent methods used there become stuck in local minima. Thus, the approach considered in [35] is reliant on the ability to find initial conditions that are close to global optima. Although it is often possible to find good initial conditions -e.g., when fitting tessellations to foams with regular structures -this is not always the case, as will be seen in Section 5. In addition, even when good initial conditions can be found, the quality of the approximating tessellation can usually be significantly improved when the optimization algorithm is able to escape local minima.
Stochastic optimization methods are widely used tools for solving highdimensional optimization problems with many local minima; see [38] for an overview. Such methods have been used to solve problems related to the LAP. In [31], the CE method is used to find a solution to the LIP and, in [33], simulated annealing is used to fit 3D Laguerre tessellations to 2D image data. However, no stochastic approach has yet been proposed in order to solve the LAP for large, noisy 3D data sets.
In our approach to solving the LAP for large tomographic data sets, we use the CE method to minimize an interface-based discrepancy measure. The CE method has been widely applied to high-dimensional multi-extremal optimization problems; see, [36][37][38][39]. It has a number of advantages over other stochastic optimization methods. For example, simulated annealing (described in [38,47]) cannot be easily applied to the LAP, as it is very difficult to find an appropriate cooling schedule. In addition, the CE method can be easily parallelized.
By using the CE method to minimize an interface-based discrepancy measure, we have a method that can escape local minima and solves a problem whose dimensionality grows linearly in the number of grains. Although the discrepancy is expensive to evaluate, we reduce this cost significantly by minimizing an interfacebased measure. In addition, we further reduce the cost of calculating the discrepancy by replacing D defined in (3) with an approximation, D.

Approximating the discrepancy
The interface-based discrepancy, D, measures the distance between voxels on the boundaries between cells in the empirical data and the planes separating the generating points in the approximating tessellation. In order to calculate this discrepancy, every voxel in each grain boundary is considered. Although this is much faster to calculate than a volume-based discrepancy, it is still computationally expensive. However, very good approximations of D can be obtained by instead considering sets of test points that describe the interfaces between the grains. We find that test points obtained by fitting approximating planes to the interfaces, using the orthogonal regression approach introduced in [48], work very well.
In order to calculate the test points, we consider only the points separating the two cells being considered (i.e., we ignore points of contact between three or more grains). Thus, instead of N i,j , we consider This avoids some numerical issues that could arise when using orthogonal regression. We then approximate the boundary between the ith and jth grain using the plane that minimizes the total squared distance to all points in N i,j . Determining this plane is a least-squares problem, which we solve using singular value decomposition. The approximating plane passes through the centroid, c ∈ R 3 , of the set N i,j . We obtain a normal vector, n ∈ R 3 , to the plane by taking the right-singular vector corresponding to the smallest singular value of the matrix containing the voxel coordinates in N i,j shifted by −c. For more details on the singular value decomposition approach, see [49]. The test points, T i,j ⊂ R 3 are chosen to lie on this approximating plane. More precisely, we put a circle with radius #N i,j /4 around the centroid c with the same orientation as the plane. We then place 10 test points equidistantly on this circle. An illustration is given in Figure 2. Illustration of test points construction. First, a circle is fitted to a grain interface. The center of the circle is given by the centroid of the voxels adjacent to both grains, its 3D orientation is obtained from the orthogonal regression plane. The radius is chosen with respect to the interface area. Finally, 10 test points are placed equidistantly on this circle.
In cases where the estimation of the plane is not 'stable' (i.e., the number of voxels in N i,j is too low to determine the location and orientation), we simply use the centroid c for all 10 test points. The criterion for 'stable' is simple: the smallest singular value must be smaller than half the second largest one. This way we can be confident that the third and shortest axis (which is perpendicular to the plane) is correctly identified.
Given the test points for the boundary between adjacent grains i and j, we define the approximate discrepancy at the boundary by where d(y, P j,k ) is given in (2) and P i,j is the plane separating the generating points (x i , r i ) and (x j , r j ). The total approximate discrepancy is then given by Note that, here, we normalize the discrepancy. This is done in order to make the cost function easier to interpret: it can be thought of as the average squared distance of test points from their corresponding separating plane.
4. Solving the LAP using the CE method

The CE method
The cross-entropy (CE) method is a stochastic optimization method that is able to solve many difficult optimization problems, including combinatorial optimization problems and continuous optimization problems with many local minima; see [36,39]. The fundamental idea of the method is to describe the location of the global minimum of an m-dimensional cost function, S(·), in terms of a degener-ate m-dimensional probability distribution. That is, a random variable with this distribution takes only a single value, namely x * = arg min x∈R m S(x). If there is more than one global minimum, each minimum will have a corresponding degenerate distribution. The CE method works by 'learning' one of these distributions.
In the continuous setting, this is done as follows. A parametric density, f (x; θ), is used to describe the possible locations of a global minimum. A sample of size M is drawn from this density and ordered by the corresponding values of S. The ρM (where · is the ceiling function) members of the sample with the lowest values of S are identified as the 'elite' sample, where ρ ∈ (0, 1). The elite sample is then used to update the parameter vector, θ. The updating is done by choosing θ to minimize the cross-entropy distance (also called the Kullback-Leibler divergence) between f (x; θ) and the targeted degenerate probability distribution. The process is continued until a stopping condition is met, e.g., the probability distribution is nearly degenerate. Thus, the CE method uses information about good choices of arguments to find better choices. The CE method and its convergence properties are discussed extensively in [36][37][38][39]. Usually, the parametric density, f (x; θ), is chosen to be a product of normal densities. That is, we choose where ϕ(·; µ, σ) is the density of a normal distribution with mean µ and standard deviation σ and θ = (µ 1 , σ 1 , . . . , µ m , σ m ). Thus, each component of the argument of S(·) is drawn independently from a normal distribution. There are two main reasons for using normal distributions. First, as the standard deviation of a normal distribution goes to zero, the distribution converges to a degenerate distribution. The second reason is that the parameter vector, θ * , which minimizes the crossentropy distance, is simply given by the maximum-likelihood estimates obtained from the elite sample. In other words, we update the parameters by setting the means equal to the sample means of the elite sample and the standard deviations equal to the sample standard deviations of the elite sample.
Using this approach, we have the following general algorithm for estimating x * = arg min x∈R m S(x).

1) Initialization. Identify an initial parameter vector,
m ) and set n = 0. 2) Sampling. Draw an independent sample X (1) , . . . , X (M ) from f (x; θ (n) ) and sort it so that S(X (1) ) ≤ · · · ≤ S(X (M ) ). We denote the jth component of Updating. Calculate the sample means and standard deviations of the elite sample. That is, calculate If a predetermined stopping condition is met, terminate. Otherwise, set n = n + 1 and repeat from step 2.

Parameter choice and stopping conditions
The parameters of the CE algorithm help to improve the speed of convergence and quality of solutions. For example, a good choice of θ (0) improves the convergence properties of the algorithm. The means, µ m , should be chosen as close as possible to the optimal solution. The initial standard deviations, σ should be chosen large enough that the algorithm is able to escape local minima but not so large that good configurations are quickly abandoned.
Ideally, the size of the sample in each step, M , should be quite large. However, there are often memory and performance constraints that limit this size. The parameter ρ, which controls the size of the elite sample, should be chosen large enough that a representative sample of good solutions are included in the elite sample. However, if ρ is chosen too large, the algorithm will take too long to converge. In general, the bigger M is, the smaller ρ should be.
A standard stopping condition for the algorithm is that it terminates when the cost function does not decrease significantly for a given number of steps.

Variance injection and dynamic smoothing
When using the CE approach outlined above, the standard deviations of the normal densities may shrink too quickly. In this case, the CE algorithm can converge to a sub-optimal solution. In order to guard against this, we use variance injection. The basic idea is to occasionally increase the variances of the distributions so that the algorithm can easily escape local minima. Usually, variance is injected when the cost function does not decrease significantly enough over a given number of iterations. The magnitude of this increase can depend on the current value of the cost function. If variance injection is used a number of times without a significant decrease in the cost function, then the algorithm is terminated.
An alternative to variance injection is to use smoothing when updating θ. The updating step for the means at the nth iteration is then of the form and the updating step for the standard deviations is given by where α, β ∈ (0, 1] (typically [0.7, 1]). It is also possible to carry out dynamic smoothing, with both α and β taken to be functions of the number of steps. In our experience, variance injection is more effective than smoothing. For a detailed discussion of both variance injection and dynamic smoothing, see [37].

Solving the LAP
We solve the LAP by minimizing the approximate discrepancy described in Section 3.4. That is, we solve where D is given in (5). Putting this in the terminology of Section 4.1, we find the arguments that minimize the cost function D, namely the N generating points, each of which is described by three coordinates and an associated radius. We associate a normal distribution with each of the m = 4N values that need to be determined. Thus, for i = 1, . . . , N , the coordinates and radius of the ith generating point, are each described by a normal density. We denote the initial means and standard deviations by µ ri . In order to ensure that the radii are positive, we truncate the corresponding normal densities to the positive real line (this is possible without changing the updating rules, see, e.g., [36,37]).

Variance injection and stopping conditions
We apply variance injection when the cost function does not decrease significantly over a period of τ > 0 iterations. More precisely, at each iteration of the algorithm, we record D (n) min , the minimal value of the approximate discrepancies calculated from the sample in that step. If, at the nth step, we perform variance injection. Because many generating points may already be close to their optimal positions and sizes, variance injection is carried out locally with a magnitude controlled by a parameter κ > 0. This is done by calculating the local cost of each cell and increasing the variance of the associated densities accordingly. We calculate the average cost of the ith cell by The local cost of a cell, D i , is defined to be the maximum of its own average cost and the average cost of its adjacent cells. That is, The variance injection is performed by setting If, at any stage, the benefit of variance injection becomes negligible, we stop performing it. More precisely, if the current minimum cost, D (n) min divided by the minimum cost immediately prior to the last variance injection is larger than γ ∈ (0, 1), we no longer carry out variance injection.
The algorithm is terminated when where δ terminate < δ inject .

Initial configuration
Reasonably good initial tessellations can be obtained directly from the tomographic data [6]. The centroids and equivalent radii of the grains in the tomographic data are used for the coordinates and radii of the generating points. The centroids, {(c xi , c yi , c zi )} i∈I , are given by

Choice of control parameters for the CE method
Our investigations showed that a parameter choice of M = 4000, ρ = 0.05, δ inject = 0.05, δ terminate = 0.01, τ = 10, γ = 0.9 and κ = 0.25 works well for a large number of different data sets. The elite set then consists of ρM = 200 samples, which results in a sufficiently large sample of good solutions. Note that other parameter values may improve the speed of convergence (or, in the case of the parameters controlling variance injection and termination, improve the quality of the approximation). However, the basic performance of the CE algorithm is relatively robust to parameter choice. That is, the convergence behavior and quality of approximations are good for a wide range of parameters.

Edge cells
The approximations obtained using the CE algorithm have unbounded cells at the edge of the observation window. This is because the corresponding grains in the tomographic data are not delimited by other grains. The unbounded cells are then intersected with the observation window, W . Note that, in most cases, edge grains are of little scientific interest and are not considered when analyzing the data. In practice, edge grains may be only partially observed or may not be representative (e.g., when stuying the dynamics of grains undergoing grain coarsening). In this paper, we explicitly ignore grains and their approximating cells if they lie at the edge of the observation window.

Experimental results
In this section, we present the results of a number of numerical experiments which we carried out in order to demonstrate the effectiveness of the CE approach. We consider three distinct data sets. Two of these data sets are produced using stochastic models. These 'artificial' data sets allow us to investigate how well our method is able to reconstruct tessellations when we are certain the underlying tessellation is, indeed, Laguerre. The first artificial data set is produced using a stochastic model that describes a polycrystalline material undergoing grain coarsening. The second artificial data set is, by design, much more pathological. It exhibits large variation in the size and shape of its grains, as well as the number of neighbors each grain has. This data set allows us to explore the effectiveness of our approach in an extreme setting. In particular, it provides an example of a setting in which the standard choice of initial conditions is far from optimal. Finally, we consider an empirical data set: tomographic data obtained from a sample of Al-5 wt% Cu. We demonstrate that our method is able to produce an excellent approximation to this data set using a Laguerre tessellation.

Artificial data
The artificial data sets are produced using two distinct stochastic models. The polycrystalline model (PCM), developed in [8], describes a polycrystalline material undergoing grain coarsening. The second data set is produced using a randomly marked Poisson process; see [50]. The resulting tessellation is known as a Poisson-Laguerre tessellation (PLT); see [43]. In the first case, there is a strong correlation between the relative positions of the generating points and their weights. In the second example, the weights are independent of the positions of the generating points. As a result, the second data set is much less regular than the first data set. Figure 3 shows 2D cross sections of the data sets, together with their generating points.

Model descriptions
In the PCM model, a random Laguerre tessellation is produced using a set of spheres that may overlap slightly. The centers of the spheres give the locations of the generating points and the radii define the weights. The spheres have hard-cores which may not overlap. This was shown to be a suitable model for polycrystalline materials in [8]. Both the density of the sphere packing and the radii of the spheres influence the sizes and shapes of the cells in the resulting tessellation. In particular, when highly dense packings are used, very "spherical" cells are generated with a narrow coordination number distribution. We generated a sample from this model in a bounded window, W = [0, 700] 3 using parameters that were fitted to a sample  of Al-5 wt% Cu annealed for 200 minutes. The fitting procedure is described in [8].
Note that the parameters we use were obtained by fitting the model to the empirical data described in Section 5.2. Having produced a sample, we removed cells that would correspond to a very small number of voxels after discretization (i.e., cells having a smaller volume than a ball with radius of 4 voxels). The final result is a tessellation consisting of roughly 2500 grains (approximately the same number as in the experimental data). The PLT model is generated by simulating a homogeneous Poisson process with some intensity λ in a bounded window W . The points are independently marked by gamma-distributed random variables with some shape parameter α and rate parameter β. Note that, unlike the PCM model, it is possible that some generating points will not produce cells in the PLT model. We generated a sample from this model in W = [0, 700] 3 with intensity λ = 7.5 × 10 −6 (so that the expected number of points is roughly 2500), where we used parameters α = 50 and β = 1.125 for the mark distribution.

Laguerre approximation
The realizations of the PCM and PLT models were both discretized on a voxel grid, resulting in two 700 × 700 × 700 images. These images are of the form described in Section 3.1. The CE method was then used to solve the LAP. The results are illustrated in Figure 4, where 2D cross-sections of the original tessellations are shown with the approximating tessellations superimposed.
Using multi-threading on a standard quad-core processor (Intel Core i5-3570K), the computing time was roughly 3 hours for the PCM data and 4 hours for the PLT data. For smaller test sets, with 500 grains each, the time required was roughly 20 to 30 minutes. The memory requirements when fitting the full data sets were minimal (especially in comparison to approaches where it is necessary to store the complete voxelized data). Basically, it is necessary to store the test points and all the generators from one iteration of the CE method. The cost of storing additional variables, such as the means and standard deviations of the densities describing the generators and the cost values of the test points are negligible. For the PCM data, the number of test points used was approximately 140 000. Approximately 160 000 points were used for the PLT data. Together with N · M ≈ 2500 · 4000 generators in one iteration, this sums up to less than 200 MB of RAM when using 32-bit floating point coordinates.   Figure 5 illustrates the convergence behavior of the CE algorithm. When approximating the PCM data, roughly 650 iterations of the algorithm were required. In the PLT case, about 900 iterations were required. It is not surprising that the CE algorithm took longer to terminate in the PLT case. As will be seen in Section 5.1.3 the initial conditions in this setting are much further from the optimal solution.

Results
In order to evaluate the quality of our approximations, we compare them to approximations obtained using the heuristic approach presented in [6] and the orthogonal regression method proposed in [48]. Note that, although the orthogonal regression approach is able to achieve quite good approximations of the cells, it does not result in a parametrized tessellation. Table 1 shows an evaluation of the approximations with respect to a volume-based discrepancy measure: the number of voxels that are  correctly labeled. Our method results in an almost perfect approximation of the PCM data. The heuristic approach considered in [6] also works very well. However, it is not able to reproduce the neighbor structure as successfully as our approach, as evidenced by rows 2 through 5 of Table 1. Our method is slightly less successful at approximating the PLT data but, in this case, considerably outperforms the heuristics of [6].
The parameters used by the heuristic approach of [6] are very close to the initial conditions of the CE algorithm. Thus, from the difference in performance, it seems that the initial conditions are quite far from the optimal parameters. In order to find better configurations, the CE algorithm has to be able to escape local minima around these initial conditions. This is made clear by the fact that increasing the number of times variance injection is used (as well as the number of iterations), we are able to further improve the results for the CE algorithm. Namely, by manually increasing the number of variance injections to 20 (instead of the 8 used in the standard approach), we obtained an approximation which correctly labeled 97.3% of the voxels (instead of 96.2%). This significant improvement indicates that local minima are present in which the CE algorithm is becoming trapped (because variance injection works by helping the algorithm escape local minima). This, in turn, implies that there are local minima near to the initial conditions of the CE algorithm that are sufficiently deep to trap it. A direct implication of these results is that algorithms that converge to nearby local minima (such as gradient-descent) may not lead to good approximations when using a standard choice of initial conditions. Furthermore, it is not clear that it is always straightforward to find good initial conditions. For example, the initial generating points are usually chosen to lie inside their cells. But, using the methods presented in [31] and the exact description of the PLT tessellation, it was not possible to find a solution to the LIP subject to the restriction that generating points were contained in their cells. If a method such as gradient-descent starts with all of the points lying inside their cells, it is unclear that it will be able to obtain a solution where some generating points lie far outside their cells. Figure 6 shows scatter plots of the original cell sizes vs. the cell sizes in the CE approximations. The quality of the CE approximation of the PCM data is immediately apparent. It is also clear that the difficulties in fitting the PLT data lie primarily in the small cells.

Experimental data
The experimental data we consider was used in [8] to develop the PCM model discussed in Section 5.1. An Al-5 wt% Cu sample (cylindrical with 8.5 mm length and 4 mm diameter) was heated to the semisolid state at a temperature of 592 • C, at which point coarsening processes were measured in situ using synchrotron X-ray tomographic imaging. For the present paper, we consider only data obtained after an annealing time of 200 minutes (which corresponds to the first annealing step in [8]). It consists of approximately 2500 grains, which is quite a large data set. A 2D cross-section of the data is shown in Figure 7. The experimental data set was segmented using the watershed transformation, resulting in a labeled image of grain   regions, which are separated by a watershed layer of one-voxel thickness. Further details regarding sample preparation, imaging and segmentation can be found in [8].

Laguerre approximation
The CE method was applied to solve the LAP for the experimental data set. A cross-section of the resulting Laguerre approximation is given in Figure 8. Fitting the approximation took approximately 70 minutes on the same Intel Core i5-3570K quad-core processor used to fit the artificial data. The number of test points extracted from the image data was roughly 150 000, which is almost the same as for the artificial data sets. The convergence behavior of the CE method is illustrated in Figure 9. Using the same parameters, only 283 iterations were necessary, which is substantially fewer than the number required for the artificial data. This is mainly due to the fact that only one variance injection was required.

Results
As with the artificial data, we compare our results with the heuristic approach proposed in [6] and the orthogonal regression approach introduced in [48]. As mentioned above, the orthogonal regression approach reconstructs the individual cells quite well but does not result in a tessellation. The results are summarized in Table 2. The CE method correctly assigns 91.4 % of the voxels. In contrast, the heuristic approach correctly assigns only 82.9 % of the voxels. The orthogonal regression yields 96.3 %. Note that the CE method significantly outperforms both the heuristic method and orthogonal regression when describing the neighborhood structure of the grains, although it does not seem possible for a tessellation containing only convex cells to accurately capture the full neighborhood structure. This is at least partially due to segmentation issues and contacts in the image data with very small areas. In contrast, the orthogonal regression approach performs surprisingly badly. It seems that the geometric properties of normal Laguerre tessellations (e.g., coinciding faces, edges and vertices) favor realistic reconstructions. We think this confirms that the Laguerre tessellation is a good choice for representing polycrystalline microstructures. It is possible, of course, that the quality of the fit could be improved using non-convex cells. For example, in [34], the linear programming method was used to obtain a generalized power diagram approximation of similar data (but with far fewer grains and voxels) resulting in a fit that correctly assigned 93.8 % of the voxels with non-convex cells.     Figure 10(b) shows that the locations of the grains are also quite accurate. Note that the main issue with fitting seems to be small cells. This is not such a problem, however, as small cells (with equivalent radii of up to 10 voxels) are subject to image segmentation error and are, thus, not too reliable in the original data.

Conclusions and outlook
In this paper, we considered the problem of approximating tomographic data by a Laguerre tessellation. We expressed this problem as an optimization problem: the generating points of the approximating tessellation need to be chosen in order to minimize the discrepancy between the tessellation and the tomographic data. We considered an interface-based discrepancy measure, instead of the volume-based discrepancies more commonly considered in the literature. This allowed us to use the CE method, a stochastic optimization method that is able to escape local minima, to fit the approximating tessellation. We then carried out numerical experiments on both artificially generated and experimentally obtained tomographic data that demonstrated the broad effectiveness of our approach.
Our method is robust and easy to implement. Thus it can be applied to fit Laguerre tessellations to almost any material with a granular or cellular structure. An obvious next step is to extend our approach to tessellations that include non-convex cells. In particular, we believe the approach and philosophy outlined in this paper can be extended to fit generalized power diagrams (see [34,51]) -generalizations of Laguerre tessellations that can include non-convex cells. This will be the subject of a forthcoming research paper.

Funding
This work was partially supported by the Australian Research Council under Grant DP140101956 and by the Deutsche Forschungsgemeinschaft under Grant KR 1658/4-1.

Supplemental material
A software package including Java code and data sets can be downloaded from https://github.com/stochastics-ulm-university/laguerre-approximation.