Laguerre tessellations and polycrystalline microstructures: A fast algorithm for generating grains of given volumes

We present a fast algorithm for generating Laguerre diagrams with cells of given volumes, which can be used for creating RVEs of polycrytalline materials for computational homogenization, or for fitting Laguerre diagrams to EBSD or XRD measurements of metals. Given a list of desired cell volumes, we solve a convex optimization problem to find a Laguerre diagram with cells of these volumes, up to any prescribed tolerance. The algorithm is built on tools from computational geometry and optimal transport theory which, as far as we are aware, are new to the microstructure modelling community. We illustrate the speed and accuracy of the algorithm by generating RVEs with user-defined volume distributions with up to 20,000 grains in 3D. We can achieve volume percentage errors of less than 1% in the order of minutes on a standard desktop PC. We also give examples of polydisperse microstructures with bands, clusters and size gradients, and of fitting a Laguerre diagram to 3D EBSD measurements of an IF steel.

In this paper we focus on the class of weighted Voronoi diagrams known as Laguerre diagrams (or power diagrams or radical Voronoi tessellations); see Definition 2.1. Laguerre diagrams have more degrees of freedom than Voronoi diagrams and so can represent more complex geometries, yet they are no more expensive to compute (Voronoi diagrams and Laguerre diagrams have the same algorithmic complexity [6]).
Goal of this paper. A limitation of Laguerre diagrams is that there is not an explicit relation between their generators and their geometric properties, such as the volumes of their cells. The goal of this paper is to develop algorithms for creating Laguerre diagrams with user-defined cell size distributions. Our motivation comes from steel modelling. We wish to generate realistic RVEs of single-and multi-phase steels for computational homogenization simulations. Unlike much of the literature on Laguerre modelling of polycrystals [1,32,33,35], our primary aim is not to fit Laguerre diagrams to EBSD or XRD data, but rather to create a tool for generating a rich family of (possibly never-observed) microstructures, which can be combined with multiscale simulations to optimize grain geometries and lead to the development of new alloys. Having said that, our algorithms are also very well suited for generating Laguerre diagrams with texture intensities that match EBSD data, as we demonstrate in Example 5.3. With these applications to steel in mind, we often refer to Laguerre cells as grains, although our results can be applied more generally to other polycrystalline metals and to foams.
Main results. Our main result is Algorithm 2, which generates 'regularized' Laguerre diagrams with grains of prescribed volumes, up to any given tolerance; see Section 3. It is an iterative method with two steps at each iteration: 1. The regularization step, which is used to avoid generating diagrams with unrealistic, highly irregular or elongated grains; 2. The optimization step, which solves a smooth, unconstrained convex optimization problem to find the weights of the Laguerre diagram so that the grains have the desired volumes. The more iterations, the more regular the Laguerre diagram.
We also provide a simplified version, Algorithm 1, which performs a single iteration with an optimization step but no regularization step. We include this for illustrative purposes (to demonstrate the effect of the regularization step of Algorithm 2) and also because it can be used to fit a Laguerre diagram to EBSD measurements of grain volumes and centroids (Example 5.3).
We present the theory underlying the algorithms in Section 2. The algorithm uses results from computational geometry and optimal transport theory [24,30], a field of mathematics that has recently gained great popularity and found applications in a wide range of areas including data science, economics, image processing, partial differential equations and statistics. We believe however that this is its first application in the steel industry.
Advantages and disadvantages of our approach. The advantages of our method are • it is fast; • it can create Laguerre diagrams with grains of exact volumes, in principle of any desired tolerance up to machine precision; • it gives some control over the spatial distribution of the grains; • it can also create periodic Laguerre diagrams, which are important for computational homogenization.
The disadvantages of our method are • it provides no direct control over the centroids of the grains or their morphology, e.g., their aspect ratio; • it is currently limited to Laguerre diagrams and so the grains cannot have curved boundaries or be non-convex.
We now discuss these advantages and disadvantages in more detail, give evidence in support of our claims, and compare our method with others in the literature.
Controlling grain volumes: Speed and accuracy. Figure 6 show that we can create Laguerre diagrams in 3D with up to 20,000 grains in around 10 minutes on a standard desktop PC (not using parallel computation), where the volumes of the grains are correct to within 1%. For 10,000 grains we require as little as 2.5 minutes (see Figure 6), although the time depends on the regularity of the microstructure; in Example 5.3 it took 6.25 minutes for 9211 grains. In our implementation of Algorithm 2 we simply used MATLAB's all-purpose fminunc optimization algorithm. With modern, customized optimal transport optimization algorithms such as [19,23] it should be possible to use our method to generate Laguerre diagrams with given volume distributions with 100,000 grains in a few minutes [23, Table 3] or even one million grains in less than an hour [23, Table 4]. The reason our method is fast is that it is equivalent to solving a sequence of convex optimization problems. We discuss this further in Section 4. We now compare this with the speed of other algorithms. It is difficult to make a direct comparison in some cases since different methods fit different geometric properties.
In [29] a stochastic optimization method (the cross-entropy method) is used to solve a nonconvex optimization problem to fit a Laguerre diagram to 3D XRD measurements of grain volumes and centroids. They report simulation times (performed using parallel computation) of 19.2 hours for 1439 grains and 122.3 hours for 8063 grains. Note that it is hard to compare their run times with ours since they are also fitting centroids; their method does not try to fit the volumes exactly like we do, but rather obtain a good fit for both volumes and centroids, and their method can produce empty Laguerre cells (grains with volume zero). While the main focus of our paper is to fit volumes only, we give an example of fitting volumes and centroids in Example 5.3, where we fit a Laguerre diagram to 3D EBSD measurements of an IF steel with 9211 grains. The run time is 6.25 minutes and the volumes are correct to within 0.56%. The centroid errors of most of the grains are comparable to those given in [29,Fig. 9], although overall our method does worse than [29] at centroid fitting, as expected. Other approaches that use stochastic optimization methods to solve non-convex optimization problems include [31,34].
Sphere-packing methods are popular for fitting Laguerre diagrams to volume distributions [12,25,28,37]. If you have n non-overlapping spheres S 1 , . . . , S n with centres x 1 , . . . , x n and radii r 1 , . . . , r n , then the Laguerre diagram with seeds x i and weights w i = r 2 i (see Definition 2.1) has the property that cell L i contains sphere S i . Therefore the volume of L i is at least the volume of S i . The idea of sphere-packing methods is that if the spheres are close-packed, then the volumes of the Laguerre cells are approximately equal to the volumes of the solid spheres. The disadvantage of this method is that it is inexact and computationally expensive since the sphere-packing problem is NP hard [18]. Nevertheless, this method provided us with inspiration for a good initial guess for the optimization simulation in Example 5.3 (see also Section 4.2.1).
In [1] a method is proposed for fitting grain measurements with generalised balanced power diagrams (GBPDs), which are a generalization of Laguerre diagrams. GBPDs are generated by triples (x i , w i , A i ) of seeds x i , weights w i , and positive definite matrices A i ; the matrices A i give some control over the aspect ratio of the generalized Laguerre cells; the case A i = I for all i corresponds to a standard Laguerre diagram. The advantage of GBPDs is that they give a high degree of control over the morphology of the grains [1,. The disadvantage is that they are hard to compute. In [1] they approximate GBPDs by voxels, and they fit discretized GBPDs to grain measurements by solving a high-dimensional linear programming problem, where the number of unknowns equals the number of grains multiplied by the number of voxels. They report that to fit a discretized GBPD to 109 grains in 3D took around 6 hours on a standard laptop (this involved solving a linear programming problem with over 77 million variables and 78 million constraints) [1,Sec. 5.3]. Again, it should be noted that it is difficult to compare their run times with ours since they are fitting grain volumes and morphology, not only grain volumes like us.
A heuristic method for approximately fitting GBPDs to measurements of grain volumes, centroids and aspect ratios was proposed by [35]. Their method does not involve solving any optimization problem at all; it includes explicit formulas for the generators (x i , w i , A i ) in terms of the data. It is reported in [35] that the method is comparable in accuracy with the optimization methods of [1,34,31] but takes a small fraction of the computation time.
No run times or volume errors are reported in [35] and so we cannot make a more precise comparison with our method.
Controlling grain geometry. The main goal of our method is to fit grain volumes fast and accurately. Unlike other methods [1,29,31,32,33,34,35] it is not specifically designed for fitting grain morphology. We now discuss to what extent we can control the geometry of Laguerre diagrams.
Our method gives some control over the spatial distribution of the grains, as shown in Figures 4,5,12, where we create microstructures with bands, clusters, and size gradients.
Several authors use grain centroids as a measure of fitting-error when fitting Laguerre diagrams to data measurements, e.g., [29,35]. We show how we can approximately fit grain centroids to 3D EBSD data in Example 5.3, although the accuracy is much lower than the volume accuracy.
In its current form our method gives no direct control over the aspect ratio of the grains. Like the sphere-packing method, Algorithm 2 tends to produce grains that are too round compared to grains typically seen in metals; see Section 4.4.
There are several ways, however, that our method could be generalised to give more control over the morphology of the grains. For example, by combining our method with multilevel Voronoi diagrams [20,39] we could maintain control over the volume of the grains while producing more realistic RVEs with non-convex and elongated grains. The idea would be to first use Algorithm 2 to create a Laguerre diagram with N 'micro-grains' of equal volume for large N . Then we would glue together the micro-grains into non-convex 'macro-grains'. By choosing which grains to glue, we would control the volume and the morphology of the macrograins. (The multilevel Voronoi method glues together two micro-grains if their generators lie in the same Voronoi cell of a 'coarser' Voronoi diagram with fewer generators.) In principle our algorithms can also be generalized very easily to produce GBPDs with grains of given volumes (up to any desired tolerance) by modifying the objective functions g and g k in Algorithms 1 and 2 (simply replace the Laguerre cells L i with generalized Laguerre cells, and replace the isotropic distances |x − x i | with anisotropic distances |x − x i | A i ). This would again allow us to control both the volumes and the aspect ratio of the grains. In practice, however, it is expensive to compute GBPDs to high accuracy; discretizing them with voxels greatly increases the cost of the algorithm. Without developing new computational geometry algorithms for the efficient computation of GBPDs, this limits the method to a small number of grains or greatly increases the run time (cf. the run time of 6 hours for 109 grains in 3D in [1]).
Since our method is currently limited to Laguerre diagrams, then the grains cannot have curved boundaries or be non-convex. Curved grain boundaries can be created using additively-weighted Voronoi diagrams (Apollonius diagrams) [6], anisotropic diagrams [3], or GBPDs [1,31,33,35], although these are all more costly to compute than Laguerre diagrams. Algorithms 1 and 2 can also be modified to produce Apollonius diagrams with grains of given volumes (in the definition of the objective functions g and g k simply replace the the Laguerre cells L i with Apollonius cells, and replace the squared distances |x − x i | 2 with non-squared distances |x − x i |) but again the implementation cost is an obstacle at the present time. We plan to explore this and the above generalizations in a future paper.
Outline of the paper. In Section 2 we recall the definition and some important properties of Laguerre diagrams. In particular Property 2.3 forms the basis of our work. Section 3 includes our main result, Algorithm 2, for generating regularized Laguerre diagrams with grains of given volumes, as well as some simple illustrative examples in 2D. Section 3 also includes Algorithm 1, which can be used for fitting a Laguerre diagram to 3D EBSD or XRD measurements of grain volumes and centroids. We discuss practical issues about how the algorithms can be implemented on a computer in Section 4. Section 5 includes some large examples (10,000+ grains) and run time tests in 3D, including examples of RVEs of Interstitial Free (IF) steels.

Laguerre diagrams
In this section we recall some facts about Laguerre diagrams.

Notation.
Let Ω ⊂ R d be the region occupied by a metal. We consider both the 2-and 3-dimensional cases (d = 2 and d = 3). For simplicity we assume that Ω is a convex polygon if d = 2 or a convex polyhedron if d = 3. In principle the algorithms below can be used for non-convex regions with curved boundaries, but they become harder to implement. In all our examples below we take Ω to be a box. If U is a subset of Ω, let |U | denote its area if d = 2 or its volume if d = 3.
Definition 2.1 ( [6,27]). Let x 1 , . . . , x n be distinct points in Ω and w 1 , . . . , w n be real numbers (not necessarily positive). The Laguerre diagram or power diagram generated by the weighted points ( We refer to the sets L i as Laguerre cells or grains.
Laguerre diagrams have the following basic properties [6,27]: • The Laguerre cells tessellate Ω, which means that n i=1 L i = Ω and cells can only intersect along their boundaries.
• If all the weights are equal, w 1 = w 2 = . . . = w n , then the Laguerre diagram is simply a Voronoi diagram.
• Adding a constant to all the weights does not affect the Laguerre diagram, i.e., the weighted points generate the same diagram for any c ∈ R.
• A generator x i need not belong to its Laguerre cell L i .
• There can be empty Laguerre cells, L j = ∅ for some j.
Now we recall two advanced properties of Laguerre diagrams, Properties 2.2 and 2.3. These are the key ingredients for generating RVEs with grains of given sizes (given areas if d = 2 or given volumes if d = 3). They are well-known in the computational geometry and optimal transport communities, but seem to be largely unknown in the microstructure modelling community. Property 2.2 states that there always exists a Laguerre diagram with grains of given sizes. Property 2.3 gives a constructive way of finding one.
The weights w i in Property 2.2 can be computed using the following result: . Then (i) The function g is concave.
. . , (x n , w n ) has cells of size m 1 , . . . , m n : Property 2.3 forms the basis of Algorithms 1 and 2. It means that if we want to generate a Laguerre diagram with grains of given sizes, then we just need to find critical points of g. Since g is concave, this is equivalent to maximizing g, or to minimizing −g, which is a smooth, unconstrained, convex optimization problem. Fast numerical methods are available for solving this [10].
Controlling the spatial distribution of grains. Property 2.3 not only allows you to control the size distribution of grains, it also gives you some control over their spatial distribution. Given positive numbers m 1 , . . . , m n with n i=1 m i = |Ω|, there are infinitely many . . , n}. This can be seen from Property 2.3; any choice of distinct points x 1 , . . . , x n gives a Laguerre diagram with grains of size m 1 , . . . , m n . In Section 4.3.1 we will show how to choose x 1 , . . . , x n to control the spatial distribution of the grains.
Connection with optimal transport theory. Properties 2.2 and 2.3 can also be stated in the language of semi-discrete optimal transport theory 1 ; see, e.g., [23], [26], [30,Sec. 6.4.2]. This connection gives us a way of finding critical points of g using fast modern methods from optimal transport theory [19,23]. We discuss this connection briefly at the end of Section 3 and will discuss it further in a future paper.

Statement of the algorithms
In this section we state the main contributions of this paper. For concreteness we state the algorithms in three dimensions, but they can also be used in two dimensions (just replace volume with area and polyhedron with polygon wherever they appear in Algorithms 1 and 2). Our main result is Algorithm 2. First however we consider a simplified version, Algorithm 1, which will help us to understand the importance of the regularization step in Algorithm 2. Algorithm 1 can also be used for data-driven modelling to fit a Laguerre diagram to EBSD or XRD measurements of grain volumes and centroids (see Example 5.3). Algorithm 1 is not new and goes back at least as far as [5]. Our role is simply to bring it to the attention of the microstructure modelling community.

Algorithm 1 Generate a Laguerre diagram with grains of given volumes
Input: A convex polyhedron Ω representing a sample of metal, a list of volumes m 1 , . . . , m n such that m i > 0 and n i=1 m i = volume(Ω), and a percentage error tolerance ε.

Method:
Initialization. Choose or randomly select n distinct points x 1 , . . . , x n in Ω. Optimization step. Use a numerical optimization method to find w = (w 1 , . . . , w n ) that maximizes the function g defined in equation (2.1). Initialize the optimization method using the initial guess w init = 0 and terminate it using the stopping criterion |∇g(w)| < 10 −2 ε min j m j . , where x = 1/185 so that the total area of all the grains equals the area of Ω. The actual areas of the grains returned by Algorithm 1 are correct to within ε = 1% error. The initialization step of Algorithm 1 was performed using the MATLAB function rand to select x 1 , . . . , x 50 (pseudo)randomly from a uniform distribution. While the grains have the correct areas to within 1%, the microstructure is somewhat irregular and unrealistic, with some highly elongated grains. This leads us to Algorithm 2, which produces more regular microstructures; compare Figures 1 and 2(i).
is the Laguerre diagram obtained in the previous iteration, which is generated by x 2. Optimization step. Use a numerical optimization method to find w (k) = w (k) 1 , . . . , w (k) n that maximizes the concave function n , w n . Initialize the optimization method using the initial guess w init = w (k−1) and terminate it using the stopping criterion ∇g k w (k) < 10 −2 ε min j m j .
n that we used for the initialization step in Example 3.1. Observe from Figure 2 how the Laguerre diagram becomes more regular as the number of iterations k increases, and how it appears to be converging. The diagram already looks quite regular after just 4 or 5 iterations and the user may be happy to take far fewer than K = 100 iterations. We discuss how to choose K in the following section.
Periodic Laguerre diagrams. Algorithms 1 and 2 can be modified to create periodic Laguerre diagrams for use as RVEs for computational homogenization (RVEs are usually taken to be periodic to avoid artificial boundary effects). To create periodic diagrams in a rectangular box Ω of side lengths l 1 , l 2 , l 3 , modify Algorithms 1 and 2 as follows. Define the periodic distance between x, y ∈ Ω by |x − y| per = min{|x − y + (il 1 , jl 2 , kl 3 )| : i, j, k ∈ Z}.
In Algorithms 1 and 2 replace the Laguerre cells L i by periodic Laguerre cellsL i , which are defined byL In Algorithm 1 replace g by Replace g k in Algorithm 2 is an analogous way.
Convergence of Algorithm 2 -centroidal Laguerre diagrams. It can be proved that Algorithm 2 converges as K → ∞, meaning that the generator locations x (k) i settle down as we take more and more iterations, like we see in Figure 2. To be more precise, there exist x 1 , . . . , x n such that lim k→∞ x (k) i = x i for all i ∈ {1, . . . , n}. We omit the proof here to avoid too many mathematical technicalities; it will appear in a forthcoming paper. We do however examine the limiting Laguerre diagram. By taking the limit k → ∞ in equation (3.1), we see that Therefore the generator x i is the centroid of its own Laguerre cell L i for all i. Such a Laguerre diagram is known as a centroidal Laguerre diagram or a centroidal power diagram, a term introduced in [11]; see also [9,38]. Centroidal Laguerre diagrams tend to be more regular than non-centroidal Laguerre diagrams, as illustrated by Figure 2(i) (centroidal) and Figure  1 (non-centroidal).  Connection with Lloyd's algorithm. If we omit the optimization step in Algorithm 2 and set the weights to be zero for all iterations, w (k) = 0 for all k, then we obtain the wellknown Lloyd's algorithm for computing centroidal Voronoi tessellations (Voronoi diagrams where each generator is the centroid of its own Voronoi cell) [13]. Therefore Algorithm 2 can be interpreted as a generalized Lloyd algorithm with capacity constraints where cell L i is constrained to have volume m i . An alternative method for generating centroidal Laguerre diagrams with capacity constraints is given in [38,Sec. 4].
Energy-decreasing property of Algorithm 2. Algorithm 2 can also be interpreted as an energy-decreasing optimization method. We describe this briefly here and postpone the details until a forthcoming paper (again, to avoid too many mathematical technicalities here). Given m 1 , . . . , m n with i m i = |Ω|, define Here the minimum is taken over all possible partitions of Ω, not just Laguerre diagrams. This is an example of an optimal transport problem. For example, in two dimensions E could represent the minimum (squared) cost of transporting the recyclable waste generated by a population uniformly distributed over a country Ω to recycling centres located at It can be shown that This is known as the optimal location problem in the economics literature [8] and the quantization problem in the discrete geometry [17], electrical engineering [16] and probability literature [15]. Thanks to the regularization step (Algorithm 2, Step 1), it can be proved that Algorithm 2 is energy-decreasing in the sense that Therefore Algorithm 2 can be viewed as a descent method for minimizing E. An alternative method for this is given in [23,Sec. 4] where, instead of updating x (k) i using our regularization step, they update it using a quasi-Newton (LBFGS) optimization step applied to E.

Computing Laguerre diagrams
One of the main expenses of Algorithms 1 and 2 is the computation of Laguerre diagrams. This happens whenever the objective function g or g k is evaluated, which could happen many times within a single optimization step. A Laguerre diagram of n generators can be computed in O(n log n) flops in 2D and O(n 2 ) flops in 3D [6, Sec. 6.2.2]. In applications n could be 10, 000 or more, and hundreds or thousands of Laguerre diagrams could be computed in a single run of either algorithm. Therefore it is important to use efficient software.
For our 2D computations we used the MATLAB function power bounded from the MAT-LAB File Exchange [40], which implements Aurenhammer's lifting method [4] and crops the diagram to a rectangular box Ω.
The function power bounded is limited to 2D, and so for our 3D computations we used (a slightly modified version of) the C++ library Voro++ [42] combined with a MEX file so that we could run Voro++ via MATLAB. In 3D we also tried the MATLAB function powerDiagramWrapper from the MATLAB File Exchange [41], combined with our own code to crop the diagram to a cuboid Ω, but we found Voro++ to be faster. Another advantage of Voro++ is that it can create periodic Laguerre diagrams.
We also used Tata Steel's own in-house Laguerre diagram code to visualise Laguerre diagrams in 2D and 3D.

Optimization methods
The other main expense of Algorithms 1 and 2 is the optimization step. For each algorithm this is a smooth, unconstrained, concave maximization problem and so is very tractable. We used the MATLAB function fminunc to minimize −g and −g k (and hence maximize g and g k ), which uses a quasi-Newton algorithm [10]. This requires an initial guess w init for the minimum point.

Choice of initial guess
For Algorithm 1 we recommend the initial guess w init = 0. For data fitting (like Example 5.3, where the seeds x i are taken from EBSD measurements), if the target grains are relatively spherical, then a better choice may be (w init ) i = m i /π in 2D or (w init ) i = (3m i /(4π)) 2/3 in 3D. In other words, (w init ) i = r 2 i where r i is the radius of a ball of area m i in 2D or volume m i in 3D. This is motivated by sphere-packing methods [12,25,28,37].
For Algorithm 2 the initial guess should depend on the iteration k. For the first iteration k = 1 we recommend w init = 0. For iterations k ≥ 2 we recommend w init = w (k−1) , the solution of the optimization step from the previous iteration. As the number of iterations increases and the points x (k) i converge, the initial guess w init = w (k−1) becomes better and better and consequently the optimization step becomes quicker and quicker. This is illustrated in Figure 3, which shows the relative run time of each iteration for Example 3.2. We see that the total runtime of the algorithm is not proportional to the number of iterations K; most of the expense is in the first few iterations.
Note that for the first iteration, k = 1, the initial guess w init = 0 does not incorporate any information about the locations x (0) i . It is possible to improve the speed of the first iteration by using a more sophisticated choice of w init , e.g., using the multilevel method of Mérigot [26] or Lévy [23], which generate a good initial guess w init by solving a sequence of smaller optimization problems with fewer grains. (For example, you can obtain a good initial guess w init for n grains by first solving a coarser problem with n/2 grains; you can obtain a good initial guess for n/2 grains by first solving a coarser problem with n/4 grains, etc.) We found that Mérigot's multilevel method [26] in 2D could halve the run time of iteration k = 1 when there are n = 10, 000 grains. It is reported that Lévy's multilevel program GEOGRAM can handle one million grains in 3D [23, Table 4]. It is also possible to obtain a better initial guess w init for iterations k ≥ 2 as follows. The Lloyd step (3.1) of Algorithm 2 could be replaced with a damped Lloyd step of the form where λ is a damping parameter between 0 and 1. The choice λ = 1 corresponds to the Lloyd step (3.1). The closer λ is to 0, the closer x , and so the better the initial guess w init = w (k−1) . Therefore the optimization step is faster for smaller λ. On the other hand, the regularization step has less effect for smaller λ, and it is necessary to increase the number of iterations K to achieve the same amount of regularization. For our purposes the full Lloyd step λ = 1 was sufficiently fast and so we did not try to optimize the choice of λ.

Choice of the tolerance
For simplicity we chose the tolerance ε of the optimization step of Algorithm 2 to be fixed at each iteration k (the optimization step terminates when ∇g k w (k) < 10 −2 ε min j m j ). The algorithm could be speeded up, however, by taking ε = ε k to depend on k. In order for Algorithm 2 to produce a Laguerre diagram with grains of given volumes up to ε% error, we only need the tolerance to be ε at the final iteration, ε K = ε. For previous iterations we could use a cruder tolerance: ε = ε K < ε K−1 < · · · < ε 2 < ε 1 . It is tempting to think that the larger the tolerance, the faster the optimization step. On the other hand, if ε k−1 is bigger than ε k , then the initial guess w init = w (k−1) at iteration k may be worse, and the optimization step at iteration k may be slower. So the tolerances ε k must be chosen carefully. The choice of fixed tolerance ε k = ε for all k is a simple, reliable option, which is why we used it.

Choice of optimization algorithm
The speed of the optimization step depends of course not only on the choice of the initial guess w init and the tolerance ε, but also on the choice of the optimization algorithm. For example, instead of using a quasi-Newton method like we did, one could use Newton's method. Newton's method tends to converge faster than quasi-Newton methods (quadratically rather than superlinearly), although it is harder to implement since it requires the second derivative of g (whereas quasi-Newton methods only require the first derivative) [10].
It can be shown (see, e.g., [9]) that where a ij is the area of the face between cell i and cell j and N i is the index set of the neighbours of cell i (that is j ∈ N i if and only if cell j and cell i share a face). The pseudoinverse of this Hessian matrix is used in a damped Newton method in [19]. The damped Newton method takes fractions of a full Newton step in order to control the error and the minimum volume of a cell (to stop cells disappearing). The method is proved to converge globally with order 1 and locally with order 2.

Initialization: Effect on the spatial distribution
In this section we discuss the initialization step of Algorithms 1 and 2.

Initialization of the seeds
The locations of the generators x 1 , . . . , x n at the termination of Algorithm 2 is a strong function of the initial choice x A further example of controlling the spatial distribution of grains can be seen in Figure  5. In this example n = 1000 grains have areas drawn from a random distribution such that   the ratio of the largest to the smallest grain size is at most 100. The Laguerre diagram in Figure 5(a) has the property that the grain sizes tend to increase form left to right. A variety of spatial distributions of grain sizes can be simulated by first distributing the seed locations appropriately. Figure 5(b) shows how a more complicated distribution can be produced.

Initialization of the weights
The choice of w (0) in the initialization step of Algorithm 2 is also important. One should choose w

Stopping criteria
In Algorithm 2 the user must specify the number of regularization steps K. As we discussed in Section 3, for large values of K the Laguerre diagram outputted by Algorithm 2 is a approximately a centroidal Laguerre diagram, which means that each seed x i . Centroidal Laguerre diagrams tend to have very regular-shaped cells, e.g., in 2D if the grains all have the same target areas, m i = 1/n, and if n and K are large, then the Laguerre diagram outputted by Algorithm 2 looks locally like a regular hexagonal tiling. Indeed for steel microstructures we found that if K is too large, then Algorithm 2 tends to produce grains that are too 'round' compared to EBSD measurements of grain aspect ratios.
Instead of fixing the number of steps K in advance, the user could terminate the algorithm whenever some measure of the maximum grain aspect ratio 2 falls below a critical threshold.
Or the user could terminate the algorithm if the distance |x | moved by the seeds from one iteration to the next falls below some threshold. The Laguerre diagram {L (k) i } tends to evolve slowly with k when k is large, as illustrated in Figure 2, and the evolution can slow down dramatically when there is a T1 topological transition (to borrow a term from foam dynamics). This topological transition involves a change of cell neighbour relations, in 2D this is via coalesence of two triple junctions of cell boundaries. So in general there is little to be gained from performing a large number of regularization steps, especially since our aim is not to produce a centroidal Laguerre diagram, but rather to produce a physically realistic microstructure.

Examples
This section includes some large examples in 3D to illustrate the capabilities of our algorithms.
Example 5.1 (Run time tests). Figure 6 gives some run times of Algorithm 2 in 3D. The black dotted line in the graph shows how the run time increases with n. We see that it grows roughly quadratically up to about n = 5000, for both the single phase case r = 1 and the dual phase case r = 5. This means that, for this range of n, the speed of the algorithm is determined by the time it takes to compute a Laguerre diagram of n cells, which is O(n 2 ) in 3D (see Section 4.1). Therefore the run time of our algorithm scales optimally in n for n < 5000 for these examples. For n > 5000 the run time grows a little faster than O(n 2 ); the cost of the optimization steps dominates the cost of computing Laguerre diagrams. Observe also from Figure 6 that it is more expensive to compute dual phase RVEs (r = 5) than single phase RVEs (r = 1).
Example 5.2 (Generating a periodic RVE of an IF steel). Figure 7 shows an example of a periodic Laguerre diagram created using Algorithm 2. The target volumes m i are taken from a 3D EBSD measurement of an IF (interstitial free) steel (EN 10130 grade DC06). There are n = 9211 grains in a box of dimensions 670µm × 80µm × 480µm. We took the initial seed locations x (0) 1 , . . . , x (0) n to be the centres of mass of the grains from the EBSD data, and we performed K = 10 regularization steps with a tolerance of ε = 0.5%. The grains in Figure  7 are coloured according to their texture (lattice orientation) by mapping the three Euler angles linearly to RGB values. The textures were taken from the EBSD data. Figure 8 shows that the volumes of all the grains are correct to within 0.5%, and most volumes are correct to within 0.05%. Example 5.3 (Fitting a Laguerre diagram to EBSD measurements). The main aim of this paper is to create Laguerre diagrams with given volume distributions for use in computational homogenization simulations. We briefly mention, however, how Algorithm 1 can be used to fit a Laguerre diagram to EBSD data of grain volumes and centroids. Figure 9 is an example of a non-periodic Laguerre diagram fitted to a 3D EBSD measurement of an IF steel (EN 10130 grade DC06) with n = 9211 grains in a box of dimensions 670µm × 80µm × 480µm (this is the same EBSD data used in the previous example). In   O(n 2 ) Figure 6: Run times in seconds of Algorithm 2 for creating single phase (monodisperse) and dual phase (polydisperse) periodic RVEs in 3D. There are n/2 grains of volume x and n/2 grains of volume rx in a cube of side length 100 with at most ε = 1% error in the volumes of the grains (r = 1 corresponds to a single phase material, r = 5 corresponds to a dual phase material, x is chosen so that the total volume of the grains equals the volume of the box). We used K = 5 regularization steps, and the initial seed locations x n were chosen randomly from a uniform distribution. The simulations were performed on an Intel Xeon E5-1620V3 (3.5GHz, 4 cores, 8 threads). The graph on the right has a log-log scale. The black dotted line is the graph of the function cn 2 , where c is a constant. It is included to illustrate how the run times grow with n.  took x 1 , . . . , x n to be the centroids of the grains from the EBSD data. The target volumes m i were also taken from the EBSD data. We used a tolerance of ε = 1%. As in the previous example, Figure 9 is coloured according to the texture of the grains. Observe that Figure 9 is less regular than Figure 7, which is due to the regularization steps of Algorithm 2. Figure 10 shows that the volumes of all the grains are correct to within 0.56%, and most volumes are correct to within 0.1%. Figure 11 shows the errors in the centroids. As expected, these are higher than the volume errors since Algorithm 1 does not directly try to fit the centroids. The relative error of 79% of the grain centroids is less than 1 and the relative error of 93% of the grain centroids is less than 2. The run time for this example was 376 seconds on an Intel Xeon E5-1620V3 (3.5GHz, 4 cores, 8 threads) with the initial guess (w init ) i = (3m i /(4π)) 2/3 , which is inspired by sphere-packing methods [12,25,28,37].
Example 5.4 (Generating a dual phase RVE with a banded microstructure). Figure 12 shows an example of a periodic, dual phase Laguerre diagram with a band of small grains in the centre, generated using Algorithm 2. There are n = 10, 000 grains: 8000 grains of volume x and 2000 grains of volume 20x (where x was chosen so that the total volume of the grains equals the volume of the box). We used K = 20 regularization steps and a volume tolerance of 1%. The grains are coloured according to their volume. In order to obtain the banded structure, we placed the initial seeds x (0) 1 , . . . , x (0) n at random within bands of the correct volume. We see from Figure 12 that these bands were largely preserved by the regularization steps.

Conclusions
In this paper we introduced a fast optimization method for generating Laguerre diagrams with user-defined grain size distributions. The volumes of the grains can be controlled exactly (to within any given tolerance). We also demonstrated how the spatial and texture distribution of the grains can be partially controlled. Our algorithms can create both non-periodic Laguerre diagrams, e.g., for data fitting, or period Laguerre diagrams, e.g., for generating RVEs of Figure 9: A non-periodic Laguerre diagram fitted to 3D EBSD data of an IF steel using Algorithm 1 (see Example 5.3). The grains are coloured according to their texture (lattice orientation). The volume distribution has a fitting error of less than 1%. The texture intensity inherits the same fitting error.  Figure 11: A histogram of the relative error of the centroids of the grains for Example 5.3. The relative error for grain i is defined by |xi − ci|/ri where xi is the centroid of grain i from the EBSD data, ci is the centroid of the Laguerre cell Li, and ri is the radius of a sphere of volume mi, where mi is the target volume of grain i. This definition of relative error was proposed by [29]. The maximum relative error is 11.16%, which is worse than the result given in [29, Fig. 9], although for most of the grains the relative errors are comparable: 7278 of the 9211 grains have relative error less than 1, and 8596 of the 9211 grains have relative error less than 2. polycrystalline metals or solid foams for computational homogenization.
A limitation of our method is that we have no direct control over the morphology of the grains, e.g., over their centroids or aspect ratio. Also our algorithms do not incorporate any physics, e.g, how the microstructure was formed during crystallization, processing, recrystallization. Consequently we found that for certain steels the grains tend to be too round. We intend to address these issues in a future paper by extending our results to a more general class of Laguerre diagrams.