On non-Poissonian Voronoi tessellations

The Voronoi tessellation is the partition of space for a given seeds pattern and the result of the partition depends completely on the type of given pattern"random", Poisson-Voronoi tessellations (PVT), or"non-random", Non Poisson-Voronoi tessellations. In this note we shall consider properties of Voronoi tessellations with centers generated by Sobol quasi random sequences which produce a more ordered disposition of the centers with respect to the PVT case. A probability density function for volumes of these Sobol Voronoi tessellations (SVT) will be proposed and compared with results of numerical simulations. An application will be presented concerning the local structure of gas ($CO_2$) in the liquid-gas coexistence phase. Furthermore a probability distribution will be computed for the length of chords resulting from the intersections of random lines with a three-dimensional SVT. The agreement of the analytical formula with the results from a computer simulation will be also investigated. Finally a new type of Voronoi tessellation based on adjustable positions of seeds has been introduced which generalizes both PVT and SVT cases.


Introduction
Three-dimensional Voronoi tessellations produce a random partition of the space which have found applications ranging from geology, (Blower et al., 2002) and molecular biology (Poupon, 2004;Dupuis et al., 2011) to numerical computing (for a review see (Du and Wang, 2005) and references therein), and chemistry (Jedlovszky et al., 2004;Idrissi et al., 2011).
In most studies Voronoi tessellations have been considered in which the positions of the centers are randomly distributed, giving rise to the so called Poisson-Voronoi tessellations (PVT) (Okabe et al., 2000) even though examples of non Poissonian Voronoi Tessellation can be found in the literature (Heinrich and Schiile, 1995;Chiu and Quine, 2001;González and Einstein, 2011). Non uniform distributions of the centers can be of interest to model regular physical configurations; we shall consider here the properties of Voronoi tessellations whose center are generated by Sobol quasi random sequences (Sobol, 1967;Bratley, P. and Fox, B. L., 1988).
A probability density function (PDF) for volumes of these Sobol Voronoi tessellations (SVT) will be proposed and compared with the results of numerical simulations.
In section 3.1 SVT and PVT will be used in an application concerning the local structure of gas (CO 2 ) in the liquid-gas coexistence phase.
In addition, we shall consider the relations between these three-dimensional structures and their lower dimensional sections; in particular we shall study chord length distributions resulting from the random intersection of SVT with straight lines. In case of PVT probability density functions of chords can be derived rigorously (Muche and Stoyan, 1992;Muche, 2010), see also (Okabe et al., 2000); here we shall present an empirical method to calculate a PDF of chord lengths in the case of SVT.
Finally Voronoi tessellations derived by perturbating the positions of points lying on a regular lattice will be considered.

Probability density functions
In the case of one-dimensional PVT, with average linear density λ, it has been proved (Kiang, 1966) that the distribution of the lengths of the segments has PDF p(l) = 4λ 2 l exp (−2λl).
(1) At the present time there are no analytical formulae for the area's and volume's distribution and we limited ourself to explore the available conjectures. Numerical experiments have shown that for PVT an approximate solution can be obtained via a 3 parameters generalized Gamma distribution, that in case of a variable x is G(x; a, b, c) = ab and the kurtosis k with Usually, instead of v the reduced variable x = v/ v is used (Tanemura, 2005): fitting procedures carried out in (Tanemura, 2005) give best parameters values a = 1.16788, b = 4.04039, c = 4.79803.
In the next section we shall adapt the Gamma three-parameter distribution to fit reduced volumes of Voronoi cells generated with Sobol sequence. The results will be compared with those obtained by a simpler PDF, a one parameter gamma used by Kiang in his seminal work on Voronoi tessellations (Kiang, 1966), whose moments are It has shown that good approximations for volume distributions of PVT cells can be obtained by setting c = 5 (Ferenc and Néda, 2007). A detailed comparison between p(x; c) and G(x; a, b, c) can be found in (Ferenc and Néda, 2007), see also (Ferraro and Zaninetti, 2012).

Volumes statistics in SVT
Sobol sequences, like all quasi-random sequences fill the space more uniformly than uncorrelated random points and this property has been extensively used in Monte Carlo methods such as integration or simulation of transport processes (Morokoff and Caflisch, 1994). Indeed the uniformity of quasi-random sequences leads to integration errors smaller than in case of random sequences. Evidence of uniformity of a Sobol sequence compared with a random one is presented in Figures 1 (155 random seeds) and 2 (140 quasi-random seeds) that show respectively a PVT and SVT. Here for the generation of points of the Sobol sequence we have used the procedure outlined in (Press et al., 1992;Antonov and Saleev, 1979), for clarity's sake the examples are two-dimensional and just 140 centers have been used. It is apparent from the Figures that SVT exhibits a narrower distribution of areas: indeed variances are 0.27 and 0.061 for PVT and SVT respectively.
A measure of uniformity of a quasi-random sequence is the discrepancy, which is the error made when representing the volume of subsets of the unit cube by the fraction of points in the subsets: the lower of the discrepancy the higher is the uniformity. There are different ways to define discrepancy see (Morokoff and Caflisch, 1994) and references therein, and it can be shown that the discrepancy on a d dimensional cube is roughly (log n) s n −1 for a large number of points n, whereas random sequences have discrepancy of size (log log n) 1/2 n −1 .
A simple measure of uniformity can be also defined as follows: for a Sobol tessellation with n seeds of an unit area the minimum distance d min between any two seeds is first computed and next d min is divided by d L = 1 n , the distance between two seeds in the case of a regular lattice of unit area. Thus the ratio dmin dL should be larger for a quasi-random sequence of seeds compared with a random one. When n=1000 the ratio dmin dL is 0.093 in the case of 2D Sobol seeds and 0.0056 in the case of random (Poissonian) seeds.
Thus, Sobol sequences present a repulsion effect that can be also found in other quasi-random sequences such as, for example, those generated by the eigenvalues of complex random matrices, see (Le Caer and Ho, 1990).
The eigenvalues seeds can be found starting from a random N × N complex matrix. The matrix elements are given by x + iy where x and y are pseudo random real numbers taken from a normal ( Gaussian ) distribution with mean zero and standard deviation 1/ √ 2. Once obtained the complex elements we diagonalize the complex matrix using the subroutine CG from the EISPACK library. The points seeds have the x and y coordinates corresponding to the real and imaginary parts of the complex eigenvalues and an example which has variance 0.06 is reported in Figure 3.
Next we have generated a three-dimensional SVT with 10 4 cells: the his-  togram of their volumes has been fitted with the generalized gamma G(x; a, b, c), as given by (2), and the one parameter gamma p(x; c), see (7), respectively. The statistical parameters of the two fits and the sample's parameters are reported in Table 1, were the first line shows numerical values of the parameters for generalized three-parameter and the one-parameter gamma, respectively, obtained from the fit of empirical data, the next lines report values of mean, variance, skewness and the kurtosis for the two distributions and the sample, here the number of cells is 10 4 and the number of bins is 40. The goodness of the fit has been assessed by first computing the PDF of both G and p and next applying the Kolmogorov-Smirnov (K-S) test (Kolmogoroff, 1941;Smirnov, 1948;Massey, 1951).
The three-parameter PDF G fits well the simulated volumes: the maximum distance between the empirical and computed distributions functions d max = 0.01 and the result of the K-S test, implemented with the FORTRAN subroutine KSONE (Press et al., 1992), gives P KS = 0.21. Note also that its moments appear to be close with those derived from the empirical distribution.
As concerns p(x:c) the results of the K-S test are worse, d max = 0.024, P KS = 1.9 10 −5 . Figure 4 shows the histogram of volumes generated by the simulation, the number of Sobol centers is 10 4 , the number of division n=40, and the graph of the generalized gamma used to fit the data with parameter values given by Table (1). In Figure 5 we report the comparison between the empirical distribution and the distribution function (DF) of G with parameters as in Table 1.
It can be interesting to compare values of the moments obtained here with those derived for PVT, by using estimate of a, b, c given in (Tanemura, 2005), namely Volumes distribution of SVT shows a smaller variance, as result of the fact the centers are distributed in a more regular fashion than in case of PVT; furthermore the PDF is more symmetric for SVT (γ SV T < γ P V T ); finally k SV T is very close to 3, the value of the Gaussian kurtosis.
A comparison between the generalized gamma PDFs for the two cases is shown in Figure 6.

An application
In this section it has been shown that the main differences between distributions of volumes in case of SVT vs PVT are that the former have a smaller variance and are more symmetric and it has been argued that, clearly, these differences are related to the more regular distributions of centers in case of SVTs. This suggests a possible application in understanding different PDF for volumes obtained in simulations of local structure of gases, in the liquid-gas phase. In (Idrissi et al., 2010) CO 2 was considered and simulations were carried out to determine the volumes available to each molecule that were considered as the center of a Voronoi tessellation. From the simulation the empirical PDF of volumes P (V ) can be computed and results show an increase of mean volume V and standard deviation σ V as temperature rises: the former effect is due to the thermal expansion of the system (Idrissi et al., 2010), whereas the increase of σ V points to more disordered distribution of the centers at higher temperatures and to increasing volume fluctuations. To compare our results with the data in (Idrissi et al., 2010), see Table 2, we have computed the standard deviation for the reduced volumes x = V / V , namely σ = σ V / V , and compared them with the standard deviations σ SV T and σ P V T of the PDFs derived in 3. As an example when T = 250 K, V = 69.7Å 3 , σ = 10.3/69.7 = 0.147, σ 2 = 0.0218 and c = 1/0.0218 = 45.79. We briefly recall that the data of the previous Table are theoretical values computed with the Voronoi polyhedra (VP) analysis based on numerical codes developed by (Jedlovszky, 1999(Jedlovszky, , 2000Tokita et al., 2004).
In the temperature range from T = 250 to T = 303 σ increases from 0.15 to 0.33 and σ SV T = 0.25 is in the middle of this range. At higher temperatures, T = 303 and T = 313 one obtains σ = 0.43 and σ = 0.41, respectively, matching closely the standard deviation derived for for PVT, namely σ P V T = 0.43. More importantly, the shape of the P DF , which is quite symmetric at the lower end of the T range, increasingly deviates symmetry as T increases and develops an exponentially decaying tail at high volume values (Idrissi et al., 2010). This is also what happens in the transition from SVT to PVT, see Figure 6. These results can be explained as follows: as T increases positions of molecules become more random and a transition takes place from a PDF relatively narrow and symmetric (like in the SVT case) to a more asymmetric PDF with larger variance corresponding to a PVT. The relation between the distributions of occupied volumes and the temperature T can be made clearer by considering the parameter c of the Kiang distribution (7). It is clear that small values of c characterize distributions with relatively large variance and skweness, whereas as c increase distributions become narrower and more simmetric. The parameter c of the Kiang function can be parameterised as function of the temperature as follows where C 1 and α 1 can be found from the data of Table 2. A numerical procedure gives C 1 = 5.7 10 25 and α 1 = −10.

Chords length distribution
In many experimental conditions it is not possible to directly observe the threedimensional cells forming a tessellation, just their linear sections: thus it is of interest to study the relationships between the geometric properties of threedimensional structures and their chords (Ruan et al., 1988;Okabe et al., 2000;Stoyan et al., 2011). One can distinguish three main ways to generate chords (Coleman, 1969), (Kellerer, 1984): isotropic uniform randomness results when the body is exposed to an uniform isotropic flow of infinite straight lines; weighted randomness occurs when a uniformly distributed random point is chosen and is traversed by a straight line with uniform random direction; two-point randomness is obtained when a straight line traverses two random points that are independently and uniformly distributed. The first case, which will be considered here since more relevant for practical applications (Kellerer, 1984); relations among PDFs of chords generated by different methods can be found in (Kellerer, 1984).
In order to obtain formulas for the distributions of the chords generated by the intersections of lines with SVTs some simplifications are needed: here the polyhedrons forming the cells will be approximated by spheres and the oneparameter distribution p as given by Equation (7), will be used to fit the cells volumes distribution. With these assumptions a formula for the PDF of chords length can been obtained, a simple iterative procedure will then be used to adapt this distribution to the simulated data, thus correcting the errors resulting from the approximations. Let p x be the probability density function for the reduced cell volumes, then the PDF p y for the lengths of the diameters is given by and the probability density function g of chords length l is see, for instance, (Ruan et al., 1988), (Watson, 1971). In the present case then, with PDF of volumes given by the Kiang function, Eq. (7), with c = 16, the distribution of diameters is p y (y) = 1 2 16 16 π 16 6 15 Γ(16) The use of the generalized gamma function, Eq.
(2), for the PDF in volumes means conversely that the integrals which follow can be done only in a numerical way. The probability density g SV T can be found by making use of Eq. (11) with p y given by (12), the result is where the coefficients are large numbers, whose values are reported in the Appendix. The previous formula and the following ones are not an exact analytical result but results from the approximation of the volume of the Voronoi's polyhedrons by spheres. In order to check the validity of this PDF we inserted in a box 50000 seeds which produce a network of irregular faces belonging to the Voronoi's polyhedra. We selected 120 random lines which will intercept the network of the irregular faces: the chord's length is evaluated as the distance between a face and the following one on the considered line. A typical run processes a total number of ≈ 3200 chords. The corresponding histogram is shown in Figure 7: note that here the results has been rescaled so that the average value of chord length is equal to 1. It is apparent that the results of the simulations do not agree with the PDF given by Eq. (13): it is enough to note that g SV T (0) = 0, is in contrast with the histogram of Figure 7. In order to overcome this problem a new variable z has been defined by a shift of l: z = l − a, so that g 1,SV T (z) = g SV T (z + a). The shift parameter a should not be confused with the parameter of the generalized gamma PDF (compare Eq. (2)). An explanation for this shift is given by the fact that the length of the chord which touches only in one point the sphere is zero conversely a chord which lies on a irregular face of the Voronoi's polyhedron has a finite length. This means that we have more short lengths in the simulation of the chords with the real polyhedrons in respect to the length's of theoretical intersections with the spheres: i.e. the PDF of having short chords is finite rather than zero. Next, to obtain a reduced variable, a scale change has been applied resulting in u = bz . The scale parameter b used here is different from the parameter  (2)). In conclusion, following translation and scale change, the final PDF is now where C is a normalizing constant, which has value 1.4717. Numerical values of a and b have been obtained by an interactive procedure that at each step computes the DF the agreement between the calculated and simulated distribution function is then been verified by the K-S test. The procedure is halted when d max , the maximum distance between the distribution functions, reaches a minimum (that is when the significance level P KS is maximum) and the corresponding pair a, b is selected; Table 4 reports the adopted values. Calculated and empirical DFs of chords length are shown in Figure 8, the K-S test gives d max = 0.0165, P KS = 0.251. The moments of g f,SV T are presented in Table 3 with parameters as in Table 4. The results can also be presented as a PDF, see Figure 9, parameters as in Table 4; the reduced χ 2 is 1.09. The shifted PDF, g f,SV T , for SVT chords can be reported as a Taylor expansion around x=0 , when the average value is one   It is clear that at x = 0 this PDF takes a finite value.
We now apply translation and scale change and Table4 reports the parameters adopted. Figure 10 reports a comparison of the previous result, integration of PDF (18)and parameters as in Table 4, as a dashed line with the tabulated result as deduced from Table 5.7.4 in (Okabe et al., 2000) when the average value of both PDFs is one. Table 5 reports the two numerical sequences.    Calculated and empirical DFs of chords length in the PVT case are shown in Figure 11 with parameters as in Table 4; the K-S test gives d max = 0.0362.
For comparison purposes a plot of PDF g f (u; a, b), which represents the SVT case with parameters as in Table 4, is shown in Figure 12 together with the numerical PDF for PVT derived from the numerical DF reported in Table  5.7.4 of (Okabe et al., 2000); in both cases the mean chord length is equal to 1. The shifted PDF, g f,P V T , for PVT chords is reported as a Taylor expansion around x = 0, when the average value is one

Adjustable seeds
An occasional reader may question if a scenario of gradual transition from PVT to SVT can be outlined. In order to have more flexible seeds we introduce the adjustable Cartesian grid (ACG) which can be computed both in 2D and 3D. The algorithm is now outlined: 1. The process starts inserting the seeds on a 2D/3D regular Cartesian grid with equal distance δ between one point and the following one
3. A random direction is chosen in 2D/3D and the two/three Cartesian coordinates of the generated radius are evaluated. These two/three small Cartesian components are added to the regular 2D/3D grid which represent the seeds. In order to have small corrections we express s in δ units. The parameter s is a good "disorder parameter" for the generated configurations. At s = 0 we will have the seeds disposed on a perfect lattice with all the volumes of the irregular polyhedra equal, increasing s we will reach before c=16 in the PDF of volumes (SVD) and subsequently c=5 (PVT). Figure 13 reports an example of 2D tessellation from ACG which areas have variance 0.047, the same value of the 2D sobol seeds. Figure 13: An example of 2D tessellation generated by 127 ACG seeds.

Applications of the adjustable seeds
We present here two applications of the adjustable seeds. These adjustable 3D seeds can be calibrated in order to have c = 16 for the Kiang distribution in volumes, see PDF (7). To each value of s corresponds one value of c that can be obtained from the relation c = 1/σ 2 . For instance when s = 0.316, c = 32.82 and at when s = 1, c = 15.87; in particular c = 16 is obtained for s = 0.64. The first application refers again to the local structure of gas (CO 2 ) in the liquid-gas coexistence phase (Idrissi et al., 2010). The parameter c of the Kiang function for the PDF in volumes can be parameterised as a function of the parameter s as follows where C 2 and α 2 can be found from a simulation. A numerical procedure gives C 2 = 24.39 and α 2 = −0.44. On equalizing the two equations (9) and (23) we obtain the following relationships between temperature, T , and regulating parameter s The previous relationship allows to find the theoretical standard deviation of the Voronoi polyhedra volumes as function of the temperature. We first fix the temperature as given by the values in Table 2 and the relationship (25) allows to find s. Given s we obtain c from eqn.(23) and by the fact that for the Kiang's function σ 2 = 1 c we easily obtain σ for a for the normalized variable. The standard deviation for the non-normalized variable is σ V = V σ. Table  6 reports the mean values and standard deviations of the Voronoi polyhedra, V , as computed in (Idrissi et al., 2010) (first line) and the procedure for the calculation of standard deviations presented here (second line). In both lines V are the values as given in (Idrissi et al., 2010) and are presented here to make the comparison easier. The third and fourth line present the values of c of the Kiang function and the regulating parameter s of the ACG seeds. Figure 14 reports the chemical standard deviation as given by Chemistry and the theoretical standard deviation as given by ACG as function of the temperature; it is clear that there is a good agreement between the calculated σ V and the experimental standard deviation. Next we consider the chords distribution of ACG tesselations, and to this end we adapt the PDF of Sobol's chords as given by (15). Calculated and empirical distribution functions of chords length for ACG with parameters as in Figure 14: Standard deviation as given by Chemistry (empty stars) and theoretical standard deviation as given by variables volumes in ACG (full points) as function of the temperature. Table 4 are shown in Figure 15; the K-S test gives d max = 0.028, P KS = 0.0051.

Conclusion
In this paper new types of three-dimensional Voronoi tessellations have been presented whose centers are not sampled from an uniform distribution, as in Poisson Voronoi (PVT) case, but rather are derived from more regular sequences. First Sobol-Voronoi Tesselations (SVT) have been considered in which cells forming the partition have as centers points generated by a Sobol quasi-random sequence. To make the notion of regularity of seeds configurations more precise a measure of uniformity has been presented. In analogy with the case of PVT we have used a generalized gamma distribution, denoted with G SV T , to fit the volumes obtained by numerical simulations. One should expect that the more regular configuration of centers in the SVT is reflected in the distributions of cell volumes and this is indeed the case: comparisons between the PDFs show that G SV T have smaller variance and are more symmetric than the distribution G P V T of volumes in respect to PVT. It should be noted that volumes distributions of PVT cells can be fitted satisfactorily by different types of gamma distributions, as mentioned in Section 2; in contrast for volumes of SVT only the generalized three-parameters gamma provides a good fit.
As concerns applications, SVT may be relevant in modelling partitions of systems that follow a more regular distribution than the usual Poisson distribution and, what is more interesting, transitions from ordered to disordered Figure 15: Data (empty circles) and theoretical DF for ACG seeds. The theoretical DF is the integral of PDF (g f,SV T ) (continuous line). states of the system must be mirrored by a corresponding change from SVT to PVT. An example has been presented in which volumes occupied by molecules of CO 2 in liquid-gas phases undergo a transformation from regular to more random distributions as the temperature increases. Transitions from regular to disordered partitions of space can be cast in a more general setting by considering the case in which center are first situated on the nodes of a regular grid and then positions are perturbed with a gaussian noise regulated by a parameter disorder s, thus creating an Adjustable Cartesian Grid (ACG). By increasing s one can reach first SVT distributions and next PVT. In this case a transition from ordered to disordered states can be parametrized by s. Considering again to local structure of CO 2 a relation has been derived between s and T , from which the standard deviation has been computed and next compared with the experimental one, showing a good agreement. Finally statistics of chords resulting from intersections of line with elements of PVT, SVT and ACG have also been investigated. The interest of such type of statistics resides in the fact that in many experimental conditions only chords of three-dimensional cells can be determined. Results show a good agreement between the analytical formula, obtained with a semi-empirical procedure and data obtained from a simulation, despite the approximations that have been used, namely considering cells to be a sphere and using an one-parameter gamma distribution to fit cells volumes, from which chords distributions have been derived. a 0 = π 2/3 8 2/3 3 √ 3, a 1 = 63997774118278000 Γ (2/3) .