Calculate the number of fixed-size circles which can fit in an arbitrary polygon (Part II)
In Calculate the number of fixed-size circles which can fit in an arbitrary polygon (Part I), I got to the point that I could generate an ideally packed collection of circles for fitting inside a figure within a fixed grid. The expansion of the hexagon to a size in which the desired circles fit with overlapping was the first step, so now we need to move on to:
- Making the circles actually overlay an input layer,
- Finding all shifted and/or rotated circle sets which have the most number of fitting circles, and
- Attempting a "best fit" among the alternatives which meet the criteria.
But first, I'm going to divert into computing theory...
It should be obvious that unless the test condition is this:
there are likely going to be an infinite number of sets of circles which meet the "most number of circles in the figure" criteria (if there are any that can). Any time the number of answers is infinite, the mechanism to obtain those answers is going to be, umm..., slow (where "slow" means "takes infinite time"). Even finding the first example of a set of circles that meets the criteria could take a while.
This brings us to "P" and "NP". Briefly, these are, respectively, the set of all "easy" problems (those that can be solved in polynomial time), and the set of "hard" problems (which can only be solved in non-polynomial time). This barely scratches the surface of the topic, but the basic theory here is that "easy" problems, where the solution complexity can be determined proportionally as a square of the number of data elements (n) being processed (as in an ugly sort solution), are different than "hard" problems, where the complexity is such that the solution is based on an exponent, like 2 to the n-th power (there is disagreement as to whether P == NP, but I'm not going to go down that rabbit hole).
Suffice it to say, heuristics are needed to solve problems which potentially have an effectively infinite number of solutions, and that perfect can be the enemy of "good enough". In this case, we probably don't need to test every possible origin and rotation of the packed circle set, and that we might be able to just generate a sample of the potential origins and rotations, to see if a "good enough" solution emerges (and if it doesn't we can just try a bit longer).
Now, there are two ways to generate the origin and rotations, systematic and random. In the systematic approach, we divide the figure evenly (which works out to be tessellations). For example, starting with a square to partition, we choose the center point (aka "reference point", level zero) then choose quarter centroids:
And then partition those (and so on):
Great, right? Well, maybe not, because our reference hexagons look like this:
with both type I and type II error (points outside and gaps within). Doh!
So maybe tessellations of a square are not going to optimally partition the hexagon space. But do we need to optimize the hexagon at all? Maybe we need to partition the circle, inscribed within the hexagon. But how do we do that? We could use even angular distance, say at 15 degrees and eighths of radius:
but that doesn't partition the space in 2D, since the gap increases with radius. If we use equidistant arc segments at each radius, we can control vertex density:
but that's sort of funky in distribution as well, and it's really hard to control the number of radial segments when trying to increase/decrease density (this works pretty well at eighths, but there was a constant in the computation, and I'd be afraid to halve or double density). I guess we could try a spiral (prorated along length, for which we can obtain a formula), but we'd still have the issue that the hexagons aren't filled.
Maybe we need to use hexagons after all. Since hexagons can be partitioned by six equilateral triangles, we could use the center as level 0, the centroids of the six triangles as level 1, then tessellate the triangles into four equilateral triangles of side s/2 (skipping the center point in each center triangle, which is already present in the level above). In Python, that works out to this simple expression (which, yes, is complicated by ternary logic):
count = 1 if levels == 0 else 6 * 4**(levels-1) + 1
Here's the points generated to level five, with triangles generated to level three:
And because it's pretty cool, I'll include the Python for allocating the points here:
import math
import arcpy
class Equilateral:
def __init__(self,p1,p2,p3,is_center=False):
self._p1 = p1
self._p2 = p2
self._p3 = p3
self.was_center = is_center
m2 = arcpy.Point((p2.X + p3.X) * 0.5,
(p2.Y + p3.Y) * 0.5)
dy = m2.Y - p1.Y
dx = m2.X - p1.X
h = math.sqrt(dx*dx + dy*dy)
cm = h * 2.0 / 3.0
theta = math.atan2(dy,dx)
self._p0 = arcpy.Point(p1.X + (math.cos(theta) * cm),
p1.Y + (math.sin(theta) * cm))
def as_wkt(self):
return ("POLYGON (({:f} {:f}, {:f} {:f}, {:f} {:f}, {:f} {:f}))".format(
self._p1.X,self._p1.Y,
self._p2.X,self._p2.Y,
self._p3.X,self._p3.Y,
self._p1.X,self._p1.Y))
def split(self):
m1 = arcpy.Point((self._p1.X + self._p2.X) * 0.5,
(self._p1.Y + self._p2.Y) * 0.5)
m2 = arcpy.Point((self._p2.X + self._p3.X) * 0.5,
(self._p2.Y + self._p3.Y) * 0.5)
m3 = arcpy.Point((self._p3.X + self._p1.X) * 0.5,
(self._p3.Y + self._p1.Y) * 0.5)
return [Equilateral(m1,m2,m3,True),
Equilateral(self._p1,m1,m3),
Equilateral(m1,self._p2,m2),
Equilateral(m3,m2,self._p3)]
def get_center(self):
return self._p0
class Hexagon:
def __init__(self,center,diameter,is_transverse=False):
radius = float(diameter) * 0.5
side = radius / math.sin(math.pi / 3)
#print("-- radius = {:f}, side = {:f}".format(radius,side))
self._p0 = center
if is_transverse:
self._p1 = arcpy.Point(center.X,
center.Y - side)
self._p2 = arcpy.Point(center.X + radius,
center.Y - side * math.sin(math.pi / 6.0))
self._p3 = arcpy.Point(center.X + radius,
center.Y + side * math.sin(math.pi / 6.0))
self._p4 = arcpy.Point(center.X,
center.Y + side)
self._p5 = arcpy.Point(center.X - radius,
center.Y + side * math.sin(math.pi / 6.0))
self._p6 = arcpy.Point(center.X - radius,
center.Y - side * math.sin(math.pi / 6.0))
else:
self._p1 = arcpy.Point(center.X + side,
center.Y)
self._p2 = arcpy.Point(center.X + side * math.cos(math.pi / 3.0),
center.Y + radius)
self._p3 = arcpy.Point(center.X + side * math.cos(2.0 * math.pi / 3.0),
center.Y + radius)
self._p4 = arcpy.Point(center.X - side,
center.Y)
self._p5 = arcpy.Point(center.X - side * math.cos(math.pi / 3.0),
center.Y - radius)
self._p6 = arcpy.Point(center.X - side * math.cos(2.0 * math.pi / 3.0),
center.Y - radius)
def as_wkt(self):
return ("POLYGON (({:f} {:f}, {:f} {:f}, {:f} {:f}, {:f} {:f}, "
"{:f} {:f}, {:f} {:f}, {:f} {:f}))".format(
self._p1.X,self._p1.Y,
self._p2.X,self._p2.Y,
self._p3.X,self._p3.Y,
self._p4.X,self._p4.Y,
self._p5.X,self._p5.Y,
self._p6.X,self._p6.Y,
self._p1.X,self._p1.Y))
def split(self):
return [Equilateral(self._p0,self._p1,self._p2,True),
Equilateral(self._p0,self._p2,self._p3,True),
Equilateral(self._p0,self._p3,self._p4,True),
Equilateral(self._p0,self._p4,self._p5,True),
Equilateral(self._p0,self._p5,self._p6,True),
Equilateral(self._p0,self._p6,self._p1,True)]
def get_center(self):
return self._p0
@classmethod
def compute_point_count(cls,levels):
return 1 if levels == 0 else 6 * 4**(levels-1) + 1
def allocate_points(self,levels):
points = [arcpy.Point(self._p0.X,self._p0.Y)]
if levels > 0:
parts = self.split()
for part in parts:
c = part.get_center()
points.append(arcpy.Point(c.X,c.Y))
for level in range(2,levels+1):
prev = parts
parts = []
for part in prev:
parts.extend(part.split())
for part in parts:
if (not part.was_center):
c = part.get_center()
points.append(arcpy.Point(c.X,c.Y))
return points
ref_hex = Hexagon(center=arcpy.Point(32.0,32.0),diameter=64)
for levels in range(0,5):
points = ref_hex.allocate_points(levels)
print("Level {:d} has {:3d} point{:s}".format(
levels,len(points),'' if len(points) == 1 else 's'))
Basically, I define two classes, Equilateral and Hexagon, and let the hexagon split itself into triangles, and let the triangles split themselves (keeping track of whether the triangle is at the center of the previous level). By returning the points as they're generated in level order, we still get good coverage if we need to mash on the Cancel button while the app is running.
I should note here that allocating those points is pretty intensive, and takes a long time for higher tessellation levels (even on a fairly recent high-end laptop):
When partitioning the rotation angle, we can skip the need to rotate a full 360 degrees (2*pi radians) and just consider 180 degrees (pi radians), because hexagons are biaxially symmetric. Here we can catch a break, because the same partitioning technique used on the square will work great on the angle, just in one dimension (with a bit of a twist): Start with zero and ninety degrees as level 0 (that is, both hexagon forms, default and transverse), then split those in half (at 45 and 135, using 180 as an implicit endpoint), and again (22.5,67.5,112.5,157.5), and again, which, to level four, looks like this:
The code for pre-calculating the rotation multiplier count doesn't even involve ternary logic:
num_rotations = 2**(levels+1)
This component is also NP, but smaller, since it's an exponent of two.
Okay, so with these partitioning levels, we can predict the number of iterations (and after some executions, generate a useful "time remaining" estimate). Here's a table showing number of iterations, and the products thereof:
The color gradations in this chart are at powers of 10 -- starting with green under ten thousand, yellow below one hundred thousand, orange under a million,..., the values in black exceed 32-bit INT_MAX (more than 2.14 billion* iterations). For obvious reasons, I'm going to cap the tessellation exponents at 9(-ish), and place a hard limit on the number of valid permutations permitted (over which a -99999 will show and a warning will be reported).
One last consideration: Systematic sampling is not the only option. We could also choose to randomly pick point placement options, using Create Random Points (Data Management). The trick here is that randomness generally requires more iterations to make sure there are no significant gaps. Here are the same number of points as generated by four hexagon tessellations (385 points) systematic first then with an equal number of random dots (overlaid, systematic crosses and random dots), then several multipliers of random densification:
What densification value would you choose? You know what, since this is all going into a tool, why not make that a parameter, too: If randomness is chosen, we can also select a densification factor (default of 3.0).
Finally, we know what the parameters are, and how many iterations it could take, and guess what? I'm going to need a Part III to go over the tool and the results thereof. (Sorry.)
- V
----
*In American usage; one thousand million in British parlance.