H Moulinec

A simple and fast algorithm for computing discrete Voronoi, Johnson-Mehl or Laguerre diagrams of points

Keywords: Voronoi diagram, Johnson-Mehl tessellation, Laguerre tessellation, micromechanics

This article presents an algorithm to compute digital images of Voronoi, Johnson-Mehl or Laguerre diagrams of a set of punctual sites, in a domain of a Euclidean space of any dimension. The principle of the algorithm is, in a first step, to investigate the voxels in balls centred around the sites, and, in a second step, to process the voxels remaining outside the balls. The optimal choice of ball radii can be determined analytically or numerically, which allows a performance of the algorithm in O(N v ln N s ), where N v is the total number of voxels of the domain and N s the number of sites of the tessellation. Periodic and non-periodic boundary conditions are considered.

A major advantage of the algorithm is its simplicity which makes it very easy to implement. This makes the algorithm suitable for creating high resolution images of microstructures containing a large number of cells, in particular when calculations using FFT-based homogenisation methods are then to be applied to the simulated materials.

Introduction

Generation of microstructures in computational homogenisation

Over the last three decades, considerable progress has been made in what might be called "computational homogenisation", as the part of computational mechanics dealing with homogenisation or micromechanics.

Until the early 1990s, the available hardware and software computer resources hardly allowed numerical simulations on realistic microstructures (whether real or simulated) of heterogeneus materials whose geometries are often very complex. Therefore, most geometries of the microstructures studied were still very simplified. To overcome these limitations, the efforts of numerical engineers in the 1990s and 2000s focused on improving existing methods such as finite elements, for example with the development of multilevel finite element method [START_REF] Feyel | A multilevel finite element method (fe2) to describe the response of highly non-linear structures using generalized continua[END_REF]), or developing new, more efficient or more adapted methods, notably the FFT-based homogenisation method introduced by [START_REF] Moulinec | A fast numerical-method for computing the linear and nonlinear mechanical-properties of composites[END_REF]. On the other hand, the progress of computer hardware, in particular the development of new parallel computers and the software tools (languages and application programming interfaces) that accompanied them, allowed numerical engineers to improve their codes, often at the cost of adapting numerical methods to this new paradigm (for example with domain decomposition methods [START_REF] Tallec | Domain decomposition methods for large linearly elliptic three-dimensional problems[END_REF], [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF], [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF])).

In parallel with research on calculation methods, studies have been carried out on the microstructures themselves. The question of the representativeness of the elementary volume studied was raised with perhaps greater acuity than before [START_REF] Gusev | Representative volume element size for elastic composites: A numerical study[END_REF], [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]). More recently, since the 2000s, the issue of the generation of artificial microstructures has received more particular attention. One of the questions is how to create an artificial microstructure that mimics a real microstructure, [START_REF] Barbe | Intergranular and intragranular behavior of polycrystalline aggregates. Part 1: F.E. model[END_REF], [START_REF] Fritzen | Periodic three-dimensional mesh generation for crystalline aggregates based on Voronoi tessellations[END_REF], [START_REF] Quey | Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing[END_REF], [START_REF] Altendorf | 3D reconstruction of a multiscale microstructure by anisotropic tesselation models[END_REF], [START_REF] Quey | Optimal polyhedral description of 3D polycrystals: Method and application to statistical and synchrotron x-ray diffraction data[END_REF], [START_REF] Van Nuland | An anisotropic Voronoi algorithm for generating polycrystalline microstructures with preferred growth directions[END_REF]). Another important issue is the efficiency of the algorithms generating the microstructures, whose performance can differ greatly from one method to another. The reader is referred to the review article by [START_REF] Bargmann | Generation of 3D representative volume elements for heterogeneous materials: A review[END_REF] which covers the topic of microstructure generation.

Context of the study

The present paper focuses on the generation of Voronoi, Laguerre and Johnson-Mehl diagrams -also called Voronoi, Johnson-Mehl and Laguerre tessellations -which are widely used for the simulation of microstructures of polycrystalline materials. More precisely, the objective being to calculate the mechanical behaviour of materials using FFT-based homogenisation methods [START_REF] Moulinec | A fast numerical-method for computing the linear and nonlinear mechanical-properties of composites[END_REF], [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]), the motivation of this study is to generate digital images of microstructures, i.e. discretisations of microstructures according to a regular grid, whose elementary data are generally called pixels in the two-dimensional case and voxels in the three-dimensional case. In the rest of this article, the term voxel is used for any dimension (2, 3 or higher). Instead of the term "digital image", which is widely used in the image processing community (Pratt (2007)), some authors prefer to speak of "raster images" as opposed to "vector images" or "vector graphics" for which the image is described by geometric shapes. The object of this study is thus referred to in the literature as "raster Voronoi tessellations" [START_REF] Lee | Raster Voronoi tessellation and its application to emergency modeling[END_REF] or "raster-based method" [START_REF] Li | A raster-based method for computing Voronoi diagrams of spatial objects using dynamic distance transformation[END_REF]), but more often as "discrete Voronoi diagrams" [START_REF] Velić | A fast robust algorithm for computing discrete Voronoi diagrams[END_REF].

A special attention is given, in the present paper, to the case of large images with a high number of grains. Indeed, it is not uncommon nowadays to apply FFT-based homogenisation methods on microstructures discretised with hundreds of millions of voxels [START_REF] Müller | Homogenization of linear elastic properties of short-fiber reinforced composites -a comparison of mean field and voxel-based methods[END_REF], [START_REF] Boittin | Numerical simulations and modeling of the effective plastic flow surface of a biporous material with pressurizedintergranular voids[END_REF], [START_REF] Chen | Analysis of the damage initiation in a sic/sic composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF], [START_REF] Vincent | Porous polycrystal plasticity modeling of neutron-irradiated austenitic stainless steels[END_REF], [START_REF] Marano | FFT-based simulations of slip and kink bands formation in 3D polycrystals: Influence of strain gradient crystal plasticity[END_REF] ) and, in the case of polycrystalline materials, composed of hundreds of thousands of grains [START_REF] Wojtacki | A micromechanical model for the secondary creep of elasto-viscoplastic porous materials with two rate-sensitivity exponents: Application to a mixed oxide fuel[END_REF]). The generation of such microstructures can lead to large computation times. For example, to generate an image of 1000 × 1000 × 1000 voxels of a Voronoi tessellation containing 1 000 000 grains, the general code Neper for polycrystal generation (Quey (2021a)) requires about 11 300 seconds on an 8-core processor Xeon Silver 4208 (2.10 GHz). Moreover, since FFT-based methods make the implicit assumption that the fields involved are periodic, we are interested, although not exclusively, in the case where the tessellations exhibit periodic conditions at the boundaries of the spatial domain studied. An additional constraint imposed by micromechanical studies is that the "sites" or "seeds" defining the tessellation are not necessarily placed on the exact positions of voxels of the image to be created. Indeed, there are some cases where the position of the sites and the discretisation grid of the image must be chosen independently, for example when several realisations of the same tessellation have to be made at different spatial resolutions.

Algorithms for creating tessellations

Simply put, a Voronoi tessellation of a given set of discrete points in a Euclidean space, hereafter referred to as sites, is a tiling of the space into regions, each of which gathers all points in the space closer to a given site than to any other. Laguerre and Jonhson-Mehl tessellations can be considered as variants of Voronoi tessellations, where the function evaluating the proximity of one point to another is different from the Euclidean distance. More precise definitions are given in section 2.

Since their introduction [START_REF] Descartes | Principia Philosophiae[END_REF], Lejeune Dirichlet (1850), [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques. deuxième mémoire. recherches sur les paralléloèdres primitifs[END_REF]), Voronoi tessellations -also called Dirichlet tessellations or Thiessen polygons -gave rise to an abundant literature. Many authors have been interested in the mathematical properties of these constructions (see [START_REF] Stoyan | Stochastic geometry and its applications[END_REF], [START_REF] Aurenhammer | Voronoi diagrams-a survey of a fundamental geometric data structure[END_REF] or [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF] for Voronoi tessellations, Møller (1992) for those of Johnson-Mehl and Aurenhammer (1991), [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF] or [START_REF] Lautensack | Random laguerre tessellations[END_REF] for Laguerre tessellations) and, in particular, in how to describe and calculate the cells of the tiling of the space that constitutes them. Numerous algorithms have been proposed, among which two families can be distinguished schematically.

1. A first category of algorithms, mainly from the robotics and computer vision communities, are interested in the calculation of the Euclidean Distance Transform (EDT), which consists in evaluating for each voxel of an image its distance to a given region of interest of the image. When these regions of interest are composed of isolated voxels, this problem reduces to that of creating an image of a Voronoi tessellation. See [START_REF] Fabbri | 2D Euclidean distance transform algorithms: A comparative survey[END_REF] for a review and comparison of existing algorithms. Several of the proposed algorithms achieve the optimal performance in O(N v ), where N v is the total number of voxels in the image. 2. Another family of methods, mainly from the computational geometry community, determines the geometric parameters of the cells (for example, the vertices, edges, facets of the polygons or polyhedra which constitute the cells of the Voronoi tessellations in the two-dimensional or three-dimensional cases, respectively). An extensive review of the existing algorithms for Voronoi diagrams, and a comparison of their efficiency can be found in [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF]. It appears that some algorithms attain the optimal time complexity O(N s ), where N s is the number of sites, but that many methods behave in O(N s ln N s ), which is the "worst-case optimal time complexity", or even worse.

However, these algorithms have a number of drawbacks for the problem addressed in this article, the statement of which was detailed in section 1.2. Although the approximations on the calculation of the distance in early EDT algorithms, imposed largely by the technical limitations of the time, have been overcome today with algorithms that achieve what is called "Exact EDT", all of the proposed algorithms, to the best of the author's knowledge, place the sites of the Voronoi cells in the exact location of a voxel, which is inadequate for the present study. For the case of interest of this study, i.e. the creation of images of tessellations, the algorithms of the second type must be completed by a step where it must be determined to which cell of the tessellation each voxel of the image belongs. So to speak, after having described the geometry of the cells of the tessellation, the algorithm must "fill" the cells with the voxels they contain. This option is rarely implemented in the existing codes but is proposed by the Neper code [START_REF] Quey | Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing[END_REF], Quey (2021a)), where it is called "raster tessellation meshing". As the whole numerical methods in question are rather delicate to implement, materials science researchers often opt either to use Neper or to implement a simple brute force algorithm whose low performance in O(N v N s ) makes it unsuitable for large simulations.

The aim of this article is to propose a simple algorithm, just a little more difficult to implement than the brute force algorithm, but with much better performance in O N v ln N s . First, in section 2 of the article, the tessellations of Voronoi, Johnson-Mehl and Laguerre are briefly introduced. In section 3, an accelerated algorithm is presented for the different types of tessellations and its performance is compared to that of Neper code, in section 4. Finally, in section 5, the algorithm is applied to some cases where the sites of the tessellations are not uniformly distributed.

Voronoi, Johnson-Mehl and Laguerre tessellations

Voronoi tessellations

A Voronoi tessellation, also called Voronoi diagram, associated with a given set S = {x s } s=1,...,N s of N s points, called "sites", "seeds" or "generators", of a given domain Ω of a Euclidean space, is a division of the domain into N s zones, called "cells", each cell s including all the points of the domain that are closer to the site s of position x s than to any other site in the set. Thus, the cell C s associated to the site s of the set S is a subdivision of domain Ω defined by

C s = {x ∈ Ω | d(x, x s ) ≤ d(x, x s ′ ), ∀s ′ ∈ S, s ′ s} , (1)
where d(x, y) denotes the distance between the two points x and y. Throughout the article, multidimensional quantities, such as vectors x and x s in relation 1, are noted in bold type while scalar data are in normal type. The distance d(, ) used above can be the usual Euclidean distance, but in the case of periodic conditions on the boundaries of domain Ω, the Euclidean distance can advantageously be replaced by what is called hereafter the L-periodic distance.

The reader is referred to Appendix A for the definition (A.6) and calculation (A.7) of the L-periodic distance. In the rest of the article, the distance d(, ) denotes either the Euclidean distance or the L-periodic distance, depending on the boundary conditions considered. In practical applications, the Euclidean space considered is most often two-or three-dimensional, but the definition of the Voronoi tessellation given above is valid for any finite dimension.

An example of a two-dimensional Voronoi tessellation, generated from 50 sites randomly arranged in the area, is shown in Figure 1. Periodic conditions on the borders of the domain have been assumed.

Johnson-Mehl tessellations

The model of the Johnson-Mehl tessellation, also called Johnson-Mehl-Avrami-Kolmogorov tessellation, has been developed to account for crystallisation in a material [START_REF] Johnson | Reaction kinetics in processes of nucleation and growth[END_REF]). To put it simply, a crystal grows from a site that appears in the volume at a given location and at a given time, until its boundaries reach those of other crystals. The tessellation is the resulting division of the space, when all the crystals have finished growing and have therefore completely covered the domain of space Ω considered. Thus, considering the set of sites given by where x s is the position of site s and t s is its birth time (i.e. the time at which the site appears in Ω and starts to grow), a given point x of the space would be reached by the growth of the crystal generated by site s -in the absence of interactions with any other site -at time t which is

S = {(x s , t s )} s=1,2,...,N s , (2)
t = d(x, x s )/G + t s , (3)
where G is the radial growth of the crystals, supposed to be constant and to be identical for all the crystals. Thus, the cell associated with a given site s is defined by

C s = {x ∈ Ω | d(x, x s )/G + t s ≤ d(x, x s ′ )/G + t s ′ , ∀s ′ ∈ S, s ′ s} (4)
or, equivalently,

C s = {x ∈ Ω | d(x, x s ) + Gt s ≤ d(x, x s ′ ) + Gt s ′ , ∀s ′ ∈ S, s ′ s} , (5)
G being supposed to be positive.

Laguerre tessellations

A Laguerre tessellation, also called power diagram, is generally presented in the literature as a partition of the space in cells defined by a set of spheres s of centres x s and radii r s . A given point x belongs to the cell C s associated with the site -or spheres which minimises the so-called "power distance" p d (x, s) defined by

p d (x, s) = d 2 (x, x s ) -r 2 s . (6)
Hence, the cell associated to the site s is defined by

C s = {x ∈ Ω | p d (x, s) ≤ p d (x, s ′ ), ∀s ′ ∈ S, s ′ s} (7)
Alternatively, from a physical point of view and as detailed in Appendix B, a Laguerre tessellation can be considered as the resulting partition of a crystallisation process starting from a set of sites S = {(x s , t s )} s=1,2,...,N s , where x s is the location of the site s and t s is its birth time, with a constant growth G of the square of the radii of the crystals. The principle remains similar to that of Johnson-Mehl tessellation: crystal s reaches a given point x at time t given by

t = d 2 (x, x s )/G + t s , (8)
and finally the cell C s associated to a given site s is defined by

C s = {x ∈ Ω | d 2 (x, x s ) + Gt s ≤ d 2 (x, x s ′ ) + Gt s ′ , ∀s ′ ∈ S, s ′ s} . (9)
This notation has the advantage of making the definitions of Johnson-Mehl and Laguerre tessellations similar.

Remark. In the particular case where all sites appear at the same time, or in the case G is zero, when using definitions 5 and 9, Johnson-Mehl and Laguerre tessellations reduce to a Voronoi tessellation.

Tessellation algorithms

In this section, algorithms are presented that solve the problem studied, the terms of which are briefly recalled here.

• The algorithms must create a digital image of a Voronoi, Johnson-Mehl or Laguerre tessellation of a set of punctual sites.

• The positions of the sites are arbitrary and do not necessarily coincide with voxel centres.

• The distance used to estimate the proximity of a voxel to a site must be exact.

• The boundary conditions of the image created can be periodic or non-periodic.

• Although most applications of Voronoi tessellation relate to two-or three-dimensional spaces, the algorithms should apply to spaces of any dimension.

• The domain Ω, in which the tessellation is created, is a d-dimensional rectangular parallelepiped, i.e. a volume of the Euclidean space of dimension d, defined by Ω = [0, L 1 ) × ... × [0, L d ) in a given orthonormal basis.

In addition, the arrangement of the sites in the volume is considered as purely random, following a uniform probability distribution. In the case of the Voronoi diagram, this type of tessellations is sometimes called "Poisson-Voronoi tessellations". The computational cost of the algorithms are estimated by counting the total number of calculations of distance between points performed during the execution of the algorithms (and in some cases, the number of calculations of the square of the distance). The effect of memory access and compilation optimisation are not considered, as it largely depends on the system on which the algorithms are implemented.

A brute force algorithm

The most simple algorithm to build an image of a Voronoi tessellation consists of calculating for each voxel its distance to all the sites and thus determining the closest site. Diagram 1 schematises this algorithm. Similar algorithms can easily be devised for Johnson-Mehl or Laguerre tessellations. In the case of Voronoi and Laguerre tessellations input :

a set of N s sites of positions {x s } in the unit cell Ω = [0, L 1 ) × [0, L 2 ) × ... × [0, L d ) output:
• an image I of n = (n 1 , n 2 , ..., n d ) voxels describing the tessellation in space domain Ω

• an image D of n = (n 1 , n 2 , ..., n d ) voxels giving the distance of each point x(k) to the closest Voronoi site

forall k ∈ {1, ..., n 1 } × {1, ..., n 2 } × ... × {1, ..., n d } do D(k) = ∞ forall s = 1 to N s do if d(x(k), x s ) < D(k) then I(k) = s D(k) = d(x(k), x s ) end end end
(but not of Johnson-Mehl tessellations), a substantial gain in efficiency can be achieved by considering the square of the distance when comparing the distances between a given voxel and all the sites, as this saves the costly calculation of a square root repeated many times. Despite this improvement, the method still performs poorly, the total cost of the algorithm being in O(N v N s ).

An accelerated Voronoi tessellation algorithm

Principle of the algorithm

Starting from the observation that the cells of a Poisson-Voronoi tessellation are all fairly similar in size, the principle of the algorithm consists first in examining a given neighborhood of all the sites, and determining for each voxel in the neighborhood whether the site considered is the closest, and then, in a second step, in processing the remaining voxels with the "brute force" algorithm described above. The algorithm is detailed below.

Detail of the algorithm

We start by choosing balls of investigation, B s , centred on each site x s of S with a given radius r 0 and a corresponding volume v 0 :

B s = {x ∈ Ω| d(x, x s ) ≤ r 0 } . ( 10
)
The choice of r 0 and v 0 is discussed in section 3.2.3. The voxels inside each investigation ball are considered and their distance to the associated site is evaluated and saved. After all the balls have been examined, two cases are possible for a given voxel of domain Ω:

• It belongs to one or more balls. In other words, the distance between the voxel and all the sites is smaller than r 0 only for the sites that are the centres of these balls. Thus the voxel belongs to the cell of the nearest of these sites.

• It belongs to no balls.

To complete the tessellation, the distance between each voxel of the second case and all the sites are evaluated and the voxel is assigned to the cell of the closest site.

Diagram 2 summarises the algorithm that performs these tasks.

Remarks.

Similarly to what has already been mentioned in section 3.1, the algorithm can be accelerated by considering the square of the distance instead of the distance.

A by-product of Algorithm 2 is the image D of the EDT (see section 1.3).

Optimal choice of the size of the investigation ball

The choice of the radius r 0 of the investigation balls plays a crucial role on the efficiency of Algorithm 2. When r 0 is large, the number of voxels inside the balls and thus the number of calculations of their distance to the centres of the balls are high. Inversely, the number of voxels outside the balls is small and becomes zero when the balls cover the whole domain Ω. When r 0 is small, the situation is the opposite: the number of distance calculations in the first step of the algorithm is small and the one of the second step is high. Between these two extremes, there is an optimal choice for the size of the investigation balls.

The number n step 1 of distance calculations in the first step of Algorithm 2 is equal to the total number of voxels inside the investigation balls. Neglecting the slight fluctuations of the number of voxels from one ball to another, due to spatial discretisation, one has

n step 1 = N v s v ′ s V , ( 11
)
where V is the volume of Ω and where v ′ s is the volume of intersection of the ball of investigation B s with the domain of interest, that is

• in the non-periodic case: v ′ s is the volume of the intersection of the ball B s and the domain Ω. The balls of investigation crossing the boundaries of Ω must be truncated from their part outside Ω, resulting in a reduction of the number of voxels to be investigated; in that case, one has v ′ s < v 0 . In the case when the ball is entirely inside the domain Ω, one has v ′ s = v 0 . And in the extreme case when the ball encompasses the whole domain Ω, one has v ′ s = V.

• in the periodic case: v ′ s is the volume of the intersection of the ball B s and the domain Ω s which is the translation of parallelepiped Ω, centred on x s . As long as the ball is smaller than the sphere inscribed in Ω s (i.e. r 0 ≤ L i /2, ∀i = 1, ..., d), one has v ′ s = v 0 . When the ball is bigger than the circumscribed sphere (i.e.

r 0 ≥ 1 2 d i=1 L 2 i ), one has v ′ s = V.
And between these two cases, v ′ s takes intermediate values.

With Voronoi tessellations, the conditions ensuring that v ′ s = v 0 (small balls compared to the domain size) are satisfied in most of the cases, and, in the periodic case, relation 11 becomes

n step 1 = N v N s v 0 V . ( 12
)
Algorithm 2: Faster algorithm to create an image of a Voronoi tessellation in a Euclidean space of dimension d input :

• a set of N s sites of positions {x s } in space domain

Ω = [0, L 1 ) × [0, L 2 ) × ... × [0, L d )
output:

• an image I of n = (n 1 , n 2 , ..., n d ) voxels describing the Voronoi tessellation in space domain Ω

• an image D of n = (n 1 , n 2 , ..., n d ) voxels giving the distance of each point x(k) to the closest Voronoi site (i.e. the map of the Euclidean Distance Transform)

Initialisation: radius r 0 of the balls of investigation is calculated using relation 17

forall k ∈ {1, ..., n 1 } × {1, ..., n 2 } × ... × {1, ..., n d } do D(k) = ∞ end
Step 1:

forall s = 1, 2, ..., N s do forall k | x(k) ∈ B s (Ball of centre x s and radius r 0 ) do if d(x(k), x s ) < D(k) then I(k) = s D(k) = d(x(k), x s ) end end end
Step 2:

forall k ∈ {1, ..., n 1 } × {1, ..., n 2 } × ... × {1, ..., n d } do if D(k) = ∞ then forall s = 1, 2, ..., N s do if d(x(k), x s ) < D(k) then I(k) = s D(k) = d(x(k), x s ) end end end end
If non-periodic conditions are assumed, one has only an upper bound,

n step 1 ≤ N v N s v 0 V , (13)
which tends to be exact when the number of grains increases and thus when the ratio of spheres crossing the boundaries of Ω decreases. In the following, expression 12 is used, as a fair approximation of the algorithm behaviour.

The number n step 2 of estimation of the distances in the second step of Algorithm 2 is given by

n step 2 = N s n out ,
n out being the number of voxels at a distance greater than r 0 from all sites. The value of n out depends on the set of sites, nevertheless, its statistics can be determined. When a site s is randomly placed in the volume, the probability for a given point of the domain to be outside the ball B s of radius r 0 centred on the site is given by

p s = 1 - v ′ s V . ( 14
)
In the same manner as above, v ′ s is approximated by v ′ s = v 0 and thus

p s = 1 - v 0 V . ( 15
)
As the positions of the N s sites are supposed to be uncorrelated, the probability for a given point to remain outside the N s balls of volume v 0 centred on the sites, is

p N s = Π s p s ,
where the notation Π s means the product for all s ∈ S. Hence

p N s = 1 - v 0 V N s .
Finally, the mean number of voxels outside all the balls is given by

E(n out ) = N v 1 - v 0 V N s
and the mean number of distance calculations in step 2 is

E(n step 2 ) = N s N v 1 - v 0 V N s .
Thus, the total number of distance calculations in Algorithm 2 is estimated by

E(n step 1+ 2 ) = N s N v v 0 V + 1 - v 0 V N s . ( 16
)
The optimal choice for v 0 , i.e. the value of v 0 which minimises E(n step 1+2 ), is therefore easy to calculate:

v 0 V = 1 - 1 N s 1 Ns-1 . ( 17
)
With this choice of v 0 one has

E(n step 1+ 2 ) = N v N s 1 + 1 N s 1 Ns-1 1 N s -1 .
When the number N s of sites is high, one has

E(n step 1+ 2 ) = N v ln N s + N v + o(N v ) . (18)
Thus, the computational performance of the algorithm is in O N v ln N s , which is significantly more efficient than the one of Algorithm 1 for which the total number of distance calculations is in O(N v N s ). When different configurations are considered for which the number of sites N s vary but the resolution, defined as the mean number of voxels per Voronoi cell, remains constant, the performance of the algorithm at constant resolution can be considered as being in O N s ln N s while the performance of the brute force algorithm is in O(N 2 s ). The number of distance calculations as a function of the radius r 0 , deduced from relation 16, is plotted in Figure 2 and compared with the actual numbers of distance calculations made during an execution of the algorithm, applied to the creation of an image of 200×200×200 voxels of a Voronoi tessellation from 1000 given sites. The correspondence is significant and the minimum of the curve is well estimated by relation 17. Expression 18, approximating the behaviour of the fast algorithm for large numbers of sites, is compared in Figure 3 with the actual number of distance calculations for the case of a 1000×1000×1000 images with varying numbers of cells, again leading to an appreciable agreement. An example of realisation of microstructure with Algorithm 2 is presented in Figure 4.

An accelerated Johnson-Mehl tessellation algorithm 3.3.1. Description of the algorithm

The principle of the accelerated algorithm remains the same for Johnson-Mehl tessellations, except that the quantity used to determine the site of the cell which a given voxel belongs to, is no longer the Euclidean distance, but the time t at which the growth of the crystal centred on the site reaches the point, given by relation 3, which can be recalled here:

t = d(x, x s )/G + t s .
The algorithm for the Johnson-Mehl tessellation thus resumes to the following: at a given time t 0 , whose choice is discussed in section 3.3.2, the spherical region of the fictitious crystallisation process of each site is examined. The radius r s of the ball of site s is given by where (.) + denotes the positive part. The voxels in the ball of radius r s centred on the site x s , i.e. the voxels satisfying

r s = G(t 0 -t s ) + , (19)
d(x, x s ) ≤ r s , (20)
are examined and the time they are reached by the crystal growth is evaluated with relation 3. After this step, there are two different cases for each voxel of the image:

• it belongs to one or more balls, i.e. relation 20 is satisfied for one or more sites, thus the voxel belongs to the cell of the one of these sites which minimises d(x, x s )/G + t s ,

• it belongs to no balls. In that case, the value of d(x, x s )/G + t s is calculated for all the sites and the voxel belongs to the cell of the site which minimises it.

The algorithm is summarised in Diagram 3. The differences with the algorithm for Voronoi tessellations are that the radii of the balls around the sites vary according to the sites, and that the quantity used to evaluate the proximity of a site is different.

Optimal choice for t 0

The value of t 0 must be chosen to minimise the computational cost of the algorithm, such as v 0 in the algorithm for Voronoi tessellations. The cost of the first step of the algorithm is proportional to the number of voxels in the balls centred on the sites, i.e.

n step 1 = N v s v ′ s V , (21)
where v ′ s is defined as in the case of Voronoi tessellation: it is the volume of the intersection of the ball B s and of the domain Ω s ( respectively Ω) in the periodic case (respectively in the non-periodic case), except that the balls B s have radii r s depending on the sites while they have an identical radius r 0 in the Voronoi algorithm. When a ball is smaller than the sphere inscribed in the domain, one has v ′ s = v s , where v s is the volume of the whole ball B s ( i.e. the volume of a ball of radius r s of relation 19). When the ball encompasses the whole domain, one has v ′ s = V . The exact expression of v ′ s , taking into account the case when the ball is between these two extreme cases is uselessly complicate and can be approximated by

v ′ s = min(V, v s ) = Vmin 1, v s V . ( 22
)
It can in no way be assumed that the balls B s are all small and that v ′ s = v s . In contrast to Voronoi tessellations, a wide dispersion of size of balls can occur, and the case when some balls reach the size of the volume must be taken into account. The cost of the first step of the algorithm thus reads

n step 1 = N v s min 1, v s V . ( 23
)
The cost of the second step depends on the number of voxels which have not been assigned to a cell in step 1, or, in other words, which are located outside all the investigation balls. The probability for a point x to be outside a given ball B s is given by p s = 1v ′ s /V , which can be approximated, using relation 22, by

p s = 1 -v s /V + . ( 24
)
Algorithm 3: Faster algorithm to create an image of a Johnson-Mehl tessellation input :

• a table of N s sites of positions and times {(x s , t s )} in space domain

Ω = [0, L 1 ) × [0, L 2 ), ... × [0, L d ) and time interval [0, T ),
• parameter G: "crystal growth rate" output :

• an image I(x) of n = (n 1 , n 2 , ..., n d ) voxels describing the Johnson-Mehl tessellation in space domain Ω temporary:

• an image D of same dimension as I Initialisation: t 0 is estimated using the method described in section 3.3.2

forall k ∈ {1, ..., n 1 } × {1, ..., n 2 } × {1, ..., n 3 } do D(k) = ∞ end
Step 1:

forall s = 1, 2, ..., N s do forall k | x(k) ∈ (Ball of center x s and radius r s = G(t 0 -t s ) + ) do if d(x(k), x s )/G + t s < D(k) then I(k) = s D(k) = d(x(k), x s )/G + t s end end end
Step 2:

forall k ∈ {1, ..., n 1 } × {1, ..., n 2 } × ... × {1, ..., n d } do if D(k) = ∞ then forall s = 1, 2, ..., N s do if d(x(k), x s )/G + t s < D(k) then I(k) = s D(k) = d(x(k), x s )/G + t s end end end end
Under the assumption that the positions and size of the different balls are uncorrelated, the probability for a given point x to be outside all the balls is approximated by

Π s 1 - v s V + . ( 25
)
The number of voxels that have to be taken into account in step 2 of the algorithm can therefore be estimated by

N v Π s 1 - v s V + , ( 26
)
and the mean of the cost of step 2 can be estimated by

E(n step 2 ) = N s N v Π s 1 - v s V + . ( 27
)
Finally, the estimation of the total cost of the algorithm is given by

E(n step 1+2 ) = N v s v ′ s V + N s N v Π s 1 - v ′ s V , (28)
with v ′ s given by expression 22 . The cost of step 1 is an increasing function of t 0 which varies from 0 when t 0 = min s (t s ) to a maximum value when the radial growth of the first site to appear (i.e. arg min s (t s ) ) exceeds the size of the domain, i.e. at a time given by

t 0 = min s (t s ) + 1 2G d i=1 L 2 i ,
in the case of periodic conditions, or at a time t 0 verifying

t 0 ≤ min s (t s ) + 1 G d i=1 L 2 i ,
in the case of non-periodic conditions. The cost of step 2 is a decreasing function of t 0 that varies from N s N v when t 0 = min s (t s ), to 0 when t 0 = min s (t s )+ 1 2G d i=1 L 2 i in the periodic case, or at a time t 0 lower than min s (t s )+ 1

G d i=1 L 2 i
in the non-periodic case. Thus the optimal value for t 0 , i.e. the value of t 0 that minimises n step 1+2 given by relation 28, must be searched between these two extreme values, e.g. with a dichotomy method, for a negligible computational cost.

Figure 5 illustrates the computational cost of the algorithm in a three-dimensional example. An image of 200 × 200 × 200 voxels of a Johnson-Mehl tessellation with 10 000 sites is created. The slight discrepancy observed for high values of t 0 between the model and the data, for G = 0.1, is most likely due to the approximation made by the model on the volume of the investigation balls when they cross the domain boundaries. This is consistent with the fact that the volume of the balls are overestimated by expression 22. Nevertheless, this does not affect the estimation of the optimal t 0 which is very satisfactory.

An accelerated Laguerre tessellation algorithm

The algorithm proposed for Laguerre tessellations is almost identical to the one of Jonhnson-Mehl tessellations, except that the investigations balls, corresponding to the crystal growth at a given time t 0 , have radii given by

r s = G(t 0 -t s ) + (29)
and that the function to evaluate the proximity of a point x to a site x s , is given by relation 8. Moreover, the optimal choice of t 0 must be searched between min s (t s ) and min s (t s )

+ 1 4G d i=1 L 2 i in the periodic case, min s (t s ) + 1 G d i=1 L 2 i in the non-periodic case.
An illustration of the behaviour of the algorithm is given in Figure 6. Simarly to the Johnson-Mehl algorithm, the optimal value of t 0 is well estimated by the model, although a slight discrepancy (here for G = 0.1) between the model and the data for high values of t 0 , which is likely due to the approximation made in the model in relation 22. An illustration of the performance of the algorithm as a function of the number of sites is presented in Figure 7. Measure of the numbers of distance calculations actually carried out in a given realisation with the algorithm (empty circles). Left: growth rate G = 0.1, the optimal choice found with the algorithm is t 0 = 1.57, right: growth rate G = 10, the optimal choice found with the algorithm is t 0 = 0.051. Optimal value of t 0 estimated following the method described in section 3.3.2 (black circle). Left: growth rate G = 0.1, the optimal choice of t 0 is t 0 = 0.426 (black circle), right: growth rate G = 10, the optimal choice of t 0 is t 0 = 0.028 (black circle).

Comparison with the Neper code

The performance of the algorithm proposed in this article is compared to that of the general code Neper (version 4.2.0.) which is widely used in the mechanics community involved in numerical simulations of crystalline materials. Neper is a versatile code developed by Romain Quey [START_REF] Quey | Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing[END_REF]) which enables to generate artificial polycrystalline microstructure of Voronoi type, and to create a mesh, typically to apply a finite element calculation, or to create an image of the tessellation, for example to visualise it or to apply a calculation using an FFT-based homogenisation method. Neper, which offers many more general features than the algorithms proposed here, can be applied to the creation of rasterised images of tessellations. Schematically, in that case, Neper proceeds in two successive steps. In a first step, starting from the position of the sites of the Voronoi tessellation, the geometry of the tessellation is determined, i.e. the vertices, edges and facets of the polyhedral cells are determined and stored, using a so-called "cell-based algorithm". In a second step, the image is created by a rasterisation process in which the cell to which each voxel of the image belongs is determined [START_REF] Quey | Centre[END_REF]). The brute force algorithm described in section 3.1, the algorithm outlined in section 3.2 and the code Neper were applied to create a three-dimensional image of a Voronoi tessellation containing a number of sites varying from 100 to 2 000 000. In each case, the same set of sites was used by the three algorithms. The dimension of the images has been fixed to 500 × 500 × 500, which is a realistic size for common calculations using an FFT-based homogenisation method. Figure 8 shows the execution times of the different algorithms, all three parallelised with OpenMP, on an 8-core processor Xeon Silver 4208 (2.10 GHz). The code Neper was applied in two configurations depending on whether only the step of creating the geometrical description of the cells was carried out (empty squares on Figure 8), or whether this step was followed by the creation of an image by "rasterisation" (black squares on Figure 8). It can be seen that the time to create the image from the geometrical description of the cells remains approximately constant at around 140 seconds for the different site numbers, but that the time taken to create the geometrical description increases with the number of sites: beyond about 10 000 sites, this time alone exceeds the time it takes to execute the accelerated Algorithm 3.2.

Unsurprisingly, the brute force algorithm has low performance and the time it takes to run is proportional to the number of sites. Nevertheless , for a low number of sites (less than 1000), and in this specific case, it can compete with Neper. The accelerated code is better in all cases and, as expected, takes an execution time of the order of the logarithm of the number of sites, which increases its relative performance as the number of sites increases. It took less than 30 seconds to create an image of 2 000 000 cells.

Application to non-uniform seed distributions

As mentioned in the introductory section 3, the algorithm was designed for sites randomly arranged in the volume according to a uniform distribution. This hypothesis plays a central role in the choice of the balls of investigation B s . Nevertheless, it is legitimate to question the performance of the method in the case this assumption is not satisfied. Two families of configurations for which the distribution of sites deviates from the uniform distribution have been considered.

In the first configuration, 12 3 (1728) sites were placed at the vertices of an hexagonal close-packed (HCP) lattice.

The volume considered

Ω = [0, 1) × [0, √ 3/2) × [0, √ 2
/3), was discretised in 1000 × 866 × 816 voxels so that the whole lattice can be obtained by periodic repetition of Ω and so that the sampling steps are very close in all directions. In a second arrangement, the sites were placed at the centers of hard spheres randomly placed in Ω and occupying a volume fraction of 50%. This ensures a minimum distance between sites of twice the radius of the hard spheres and thus a more ordered arrangement than the purely random distribution. In the third arrangement, the sites were placed randomly in the volume. The time spent by Algorithm 2, on an 8-core processor Xeon Silver 4208 (2.10 GHz), to build the correponding Voronoi tessellations was respectively 20.7s, 21.3s and 22.0s, to be compared to the 902.2s needed by the brute force algorithm. Algorithm 2 performs very well in configurations 1 and 2, it behaves even slightly better than in the pure random case. This is not actually very surprising as the investigation balls of the algorithm cover very well the domain Ω in those cases with an optimal coverage with the HCP arrangement. A deeper analysis of the behaviour of the algorithm shows that in the HCP and close to HCP configurations, no voxels are processed by the step 2 of the algorithm instead of 7.5% in the pure random configuration.

In the second family of configurations, the arrangements considered mimic columnar grain growth in polycrystals. In the first configuration (Figure 9(a)), 1000 sites were placed randomly in the plane z = 0 and a Voronoi tessellation was build in a 3D domain Ω = [0, 1) × [0, 1) × [0, 1), discretised in an image of 500 × 500 × 500 voxels. This a worstcase situation as the investigation balls centered around the sites cover only a small part of the whole domain, and the number of voxels processed by the brute force method in step 2 of Algorithm 2 remains high. The time spent to build the tessellation on an 8-core processor Xeon Silver 4208 (2.10 GHz) was 95.5 seconds which is comparable (and even worse) to the 92.7 seconds required by the brute force algorithm, while the same configuration of sites used to build a 2D tessellation required only 0.015 seconds. But it can be opposed that this configuration is unrealistic, as it would be much more appropriate to perform the calculation in 2D, and, if a 3D image is to be eventually constructed, to replicate the 2D image as many times as there are voxels in the third dimension. A second configuration is proposed, more realistic as it cannot be reduced to a 2D calculation: 1000 sites are randomly placed in the sub-volume z ∈ [0, 0.1) of the whole domain Ω = [0, 1) × [0, 1) × [0, 1), the latter being discretised by 500 × 500 × 500 voxels (Figure 9(b)). The calculation time on the same computer is 95.5s to be compared to the 92.7s with the brute force algorithm and the 4.8s needed by the calculation of a tessellation in the volume [0, 1) × [0, 1) × [0, 0.1) discretised in 500 × 500 × 50 voxels. A much more favourable configuration, consisted in adding to the 1000 sites in the sub-volume z ∈ [0, 0.1), 1000 other sites randomly placed in the complementary sub-volume (z ∈ [0.1, 1)), resulting in a concentration of sites 9 times lower in the region z ∈ [0.1, 1) than in the region z ∈ [0, 0.1) (see Figure 9c). The time needed to build the Voronoi tessellation was 6.5s.

In conclusion of this short study on non-uniform distributions of sites, the proposed algorithm performs very well when the sites are well distributed throughout the volume. In the opposite case, i.e. when there are both areas with high concentrations of sites and areas without any sites, the performance of the algorithm remains similar to that of the brute force algorithm.

Conclusion

In this paper, an algorithm has been presented for creating digital images of Voronoi, Johnson-Mehl and Laguerre tessellations of punctual sites assumed to be randomly distributed in the domain considered. The principle of the algorithm is to assign the voxels belonging to balls centred around the sites of the tessellation, to the cell corresponding to the closest of these sites. The choice of the radii of these investigation balls plays a predominent role on the performance of the algorithm. An analytical development enables to find a close to optimal choice, whose validity has been verified by numerical tests.

Incidentally, a relevant definition of the distance is introduced for the case of a periodic microstructure and a simple expression is given to calculate it. Thus, building an image of a unit cell under non-periodic conditions, consists in determining, for each of its voxels, if it belongs to a ball (or disc in 2D) by comparing the Euclidean distance between the voxel and the centre of the ball to the radius of the latter. In the case of periodic conditions, one has to compare the distance between the location x of each voxel and every periodically repeated centres of the balls, to the radius. In other words, one has to evaluate where x i and y i are the components of vectors x and y in an orthonormal basis. In the particular case where the sites are inside the unit cell and that the ball diameter is smaller than the dimensions of the unit cell, it is sufficient to estimate the distances between the voxels of the image of the unit cell and the ball centres which has been repeated by one period forward and backward in every directions [START_REF] Fritzen | Periodic three-dimensional mesh generation for crystalline aggregates based on Voronoi tessellations[END_REF], [START_REF] Bargmann | Generation of 3D representative volume elements for heterogeneous materials: A review[END_REF] (each component k i of k only takes the value -1, 0 or 1). Thus, to evaluate expression A.5, one has to calculate 9 Euclidean distances in a two-dimensional problem, and 27 Euclidean distances in a three-dimensional problem. However, this approach does not work in the more general case where the points whose distance is to be calculated can be located outside the unit cell. This typically occurs in the case of a microstructure that contains inclusions whose size exceeds one of the dimensions of the unit cell (for example long fibers in a periodic microstructure). Therefore, we define the L-periodic distance between two points of coordinates x and y of the Euclidean space -which can be located outside the unit cell -where the property considered is supposed to be periodic, with periods L i , with For a given growth constant G, one can define a time t s for each site s as t s = -r 2 s /G and thus each cell C s can be explicited as

C s = {x ∈ Ω | d 2 (

Figure 1 :

 1 Figure 1: Voronoi tessellation of 50 punctual sites, with periodic conditions. The cells are plotted in arbitrary colours. The sites are marked by small black circles.

Algorithm 1 :

 1 Brute force algorithm to create an image of a Voronoi tessellation from a set of positions of sites in a Euclidean space of dimension d. k gathers the d indices of the voxels in the image. x(k) is the location of voxel k.

Figure 2 :

 2 Figure 2: Number of calculations of distances between voxels and sites per voxel and per site for different values of parameter r 0 (radius of the balls of investigation around the sites), performed during the generation of an image of 200 × 200 × 200 voxels of a Voronoi diagram of 1000 sites with Algorithm 2. Empty circles: actual number of distance calculations during the execution of the algorithm. Dashed line: model given by relation 16. Black circle: optimal radius r 0 ≃ 0.12, estimated using relation 17.

Figure 3 :

 3 Figure 3: Number of calculations of distances between voxels and sites in the fast algorithm in the case of an image of 1000 × 1000 × 1000 voxels of a Voronoi tessellation for different numbers of sites N s . Measure of the numbers of distance calculations actually performed during the execution of the algorithm (circles). Expression 18 approximating the behaviour of the algorithm for large numbers of cells (dashed line).

Figure 4 :

 4 Figure 4: Image of 1000 × 1000 × 1000 voxels of a Voronoi tessellation of 1 000 000 punctual sites, generated by Algorithm 2 in 77.2 seconds on an 8-core processor Xeon Silver 4208 (2.10 GHz). The cells are plotted in arbitrary pseudo-colours. Left: general view of the volume. Right: detail.

Figure 5 :

 5 Figure 5: Number of calculations of distances between voxels and sites made in Algorithm 3 to generate an image of 200 × 200 × 200 voxels of a Johnson-Mehl tessellation of 10000 sites. Model given by relation 28 (continuous line) with varying values of t 0 .Measure of the numbers of distance calculations actually carried out in a given realisation with the algorithm (empty circles). Left: growth rate G = 0.1, the optimal choice found with the algorithm is t 0 = 1.57, right: growth rate G = 10, the optimal choice found with the algorithm is t 0 = 0.051. Optimal value of t 0 estimated following the method described in section 3.3.2 (black circle).

Figure 6 :

 6 Figure 6: Number of calculations of distances between voxels and sites made in the algorithm to generate an image of 200 × 200 × 200 voxels of a Laguerre tessellation of 10000 sites. Model giving by relation 28 and 29 (dashed line) with varying values of t 0 . Measure of the numbers of distance calculations actually performed during the execution of the algorithm (empty circles).Left: growth rate G = 0.1, the optimal choice of t 0 is t 0 = 0.426 (black circle), right: growth rate G = 10, the optimal choice of t 0 is t 0 = 0.028 (black circle).

Figure 7 :

 7 Figure 7: Total execution times t e of the accelerated algorithm for the creation of an image of 1000 × 1000 × 1000 voxels of a Laguerre tessellation (radial growth G = 0.01) for different numbers of sites N s on an 8-core processor Xeon Silver 4208 (2.10 GHz) Dashed line: least square regression on the data ln N s , t e .

Figure 8 :

 8 Figure 8: Code performance comparison on an 8-core processor Xeon Silver 4208 (2.10 GHz). Generation of an image of 500 × 500 × 500 voxels of a Voronoi tessellation with different numbers of sites. Total execution time in seconds for the brute force algorithm (empty circles), for the accelerated algorithm (black circles), for Neper code (black squares), for Neper code without image rasterisation (empty squares).

Figure 9 :

 9 Figure 9: 3D Voronoi tessellations simulating columnar grain growth in a polycrystal. (a): 1000 sites randomly placed in the plane z = 0. The tessellation was build in the 3D volume Ω = [0, 1) × [0, 1) × [0, 1). (b): quasi-columnar growth, 1000 sites were randomly placed in the subvolume z ∈ [0, 1) of Ω. (c): 1000 sites were randomly placed in the subvolume z ∈ [0, 1) of Ω, and 1000 other sites in the subvolume z ∈ [0.1, 1). The 3 images were discretised into 500 × 500 × 500 voxels.

Figure

  Figure A.11: Left: unit cell of the two-dimensional microstructure A.10. Right: image of the same discs without periodic conditions.

  d ′ (x, x s ) = min k∈Z d d E (x, x s + d i=1 k i L i ) , (A.3)where x denotes the position of the voxel and x s the one of the site or centre, where k is a d-dimensional vector of integer components k i , and where d E denotes the usual Euclidean distance:d E (x, y) = d i=1 (x iy i ) 2 ,(A.4)

  [START_REF] Yan | Computing 2D periodic centroidal Voronoi tessellation[END_REF], and one hasd ′ (x, x s ) = min k∈{-1,0,1} d d E (x,

  x, x s ) + Gt s ≤ d 2 (x, x s ′ ) + Gt s ′ , ∀s ′ ∈ S, s ′ s} . (B.3) The time t s at which a site "appears", if one takes again the analogy with the crystallisation process, does not necessarily have to be negative if one remarks that the partition obtained does not depend on the choice of the time reference, i.e. the partition is unchanged if one considers time t ′ = tt re f instead of t, the power-distance being modified intod 2 (x, x s ) + Gt ′ s = d 2 (x, x s ) + Gt s -Gt re f (B.4) because min s d 2 (x, x s ) + Gt ′ s = min s d 2 (x, x s ) + Gt s -Gt re f . (B.5)The relation between the radius r s and time t s can then be expressed asG(t st re f ) = -r 2

s (B.6) with t re f chosen adequately.

Acknowledgments

The author thanks the "Institut de Radioprotection et de Sûreté Nucléaire", and especially Dr Pierre-Guy Vincent, who has supported the project for which the method proposed in this article has been developed.

The author also thanks Romain Quey for his interest in this study and for his pertinent advice.

Available code

A beta version of the code applying the algorithms presented in this article is freely available on the site: https://lma-software-craft.cnrs.fr/voronoi-johnson-mehl-and-laguerre-tessellations .

Appendix A. L-periodic distance

Appendix A.1. Periodicity In the Euclidean space of dimension d, a given property f defined at each point x is periodic of periods L i , with i = 1, 2, ..., d, when

where L i are given d-dimensional vectors supposed to be linearly independent of each other. The property f is completely determined by its definition in the unit cell

or in any translation of U by a vector o:

. Definition of the L-periodic distance

To introduce the definition of the L-periodic distance, the simple example of a two-dimensional microstructure containing two circular inclusions periodically repeated in a matrix is considered. An illustration of the resulting microstructure is presented in Figure A.10 where a period of the pattern is repeated three times in each direction. The

. Calculation of the L-periodic distance

In this paragraph, the expression of the L-periodic distance is developed to obtain a formulation convenient to use in practice.

It is straightforward that relation A.6 simplifies into

In the case when the vectors L i are all orthogonal to each other, they can be noted

where L i are positive scalars and e i are the unit vectors of an orthonormal basis. Hence, d L can be

where x i and y i are the components of x and y in the basis (e i ). With no real difficulties, it can be established that

where N() is the "nearest integer function" defined as

where ⌊.⌋ denotes the "floor function", that gives the greatest integer value lower or equal than its real number argument. The nearest integer function can be calculated using the function rint() of the C mathematical library. Finally, the L-periodic distance d L can be written

The numerical cost of the calculation of the L-periodic distance using expression A.7 is only slightly higher than the one of the Euclidean distance with expression A.4. In addition to being an exact method, it is very simple to implement and its numerical cost is much lower than evaluating the minimum of 27 (respectively 9) Euclidean distances in 3D (respectively 2D) problems, as proposed by [START_REF] Fritzen | Periodic three-dimensional mesh generation for crystalline aggregates based on Voronoi tessellations[END_REF] and [START_REF] Bargmann | Generation of 3D representative volume elements for heterogeneous materials: A review[END_REF].

Appendix B. Alternative definition of Laguerre tessellations

As already presented in section 2.3, a Laguerre tessellation associated with the set of sites S = {(x s , r s )} s=1,2,...,N s where x s is the position of site s and r s is an associated positive scalar, is a division of the domain Ω of the Euclidean space into N s cells defined by