circle_packing, a FORTRAN90 code which analyzes point sets in the unit square (or M-dimensional unit hypercube).
The point sets are to be measured in terms of a circle-packing test. The point of the test can best be described in analogy. Suppose that, at each point in the set, we place a (two dimensional) balloon of initial radius 0. At a given instant, all the balloons begin to expand. If any two balloons touch, they instantly stop expanding. Once all the balloons have stopped expanding, we record the radius of each balloon.
The point sets we are interested in have been generated by a variety of algorithms, and are stored as datasets, of spatial dimension 2, 7 and 16, and sample sizes 10, 100, 1,000 and 10,000. The point set types include:
The amount of area covered by the circles is compared to the area of the unit square. This measurement has a certain amount of boundary effect, and is more effective when the number of points is large, so that the boundary effects can be ignored. In order to make a meaningful number to report, we prefer, then, to restrict the balloons so that they cannot extend outside the unit square. In this way we can compare the coverage for a given set of points to that produced by a maximal sphere packing set, for which the balloons are all of the same size.
Note that, for each datapoint, the radius of the associated circle is the minimum distance to the nearest boundary of the Voronoi cell containing the point.
Based on this measure of coverage, a "good" dataset would be one for which the value was maximal. We might also be interested in the variation in the sizes of the radiuses.
In the limit, the maximum relative packing density of a 2D region with equal sized spheres is 0.9069. In 3D, a density of at least 0.74 can be achieved, and it is known that no greater than 0.7796 is possible.
(in)Efficiency Note: As currently programmed, the routine is inefficient, taking execution time that is O(N^3). This is because there are about N steps of the iteration, each of which "retires" one point by finding its maximal radius; in each of the N steps, the program recomputes the distances between every unretired point and all the N points. A more clever program would construct the Voronoi diagram, and then each point would only have to worry about its nearest neighbors.
For instance, for the 2D Halton datasets, here is the table of the number of points versus unrestricted (spheres can extend out of the unit hypercube) and restricted (they can't) scores:
| Points | Unrestricted | Restricted |
|---|---|---|
| 10 | 0.8203 | 0.3884 |
| 100 | 0.5416 | 0.4067 |
| 1000 | 0.4299 | 0.3959 |
| 10000 | 0.4214 | (unavailable) |
A collection of sample datasets, and plots of the points and the circles involved in the circle test, is available at sample_2d.
circle_packing input_filewhere
The computer code and data files made available on this web page are distributed under the MIT license
circle_packing is available in a FORTRAN90 version.
DIAPHONY, a C program which reads a file of N points in M dimensions and computes its diaphony, a measure of point dispersion.
STAR_DISCREPANCY, a C program which reads a file of N points in M dimensions, (presumed to lie in the unit hypercube) and computes bounds on the star discrepancy, a measure of dispersion, by Eric Thiemard.
TABLE_QUALITY, a FORTRAN90 code which computes the quality measures of a dataset read from a file. which can analyze a dataset that is stored in a file.
TET_MESH_QUALITY, a FORTRAN90 code which computes quality measures of a tetrahedral mesh.
TRIANGULATION_QUALITY, a FORTRAN90 code which computes quality measures of a triangulation.