Latest Contributions by VinceAngelo

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. ... View more

I've seen a number of queries recently, trying to find a tool which would allow a user to determine the maximum number of circles of a certain size which fit inside an arbitrary polygon.  There is, of course, no such COTS tool, but generating a custom tool task starts with generating a region of tightly packed circles, and there's almost a tool to do that, though it needs some tweaking... The ArcGIS GenerateTessellation tool (arcpy.management.GenerateTessellation) will generate tessellations of a number of given shapes*: *Note: Transverse Hexagon and Diamond only available in ArcGIS Pro but "Circle" isn't included because the definition of tessellation ("a covering of an infinite geometric plane without gaps or overlaps by congruent plane figures of one type or a few types" -- Merriam-Webster) requires the figures to be "without gaps or overlaps". However, if your intent is to use circles anyway, then there are clearly two base geometries one could use to start: SQUARE (which is the same as DIAMOND, just rotated 45 degrees), and HEXAGON (for which TRANSVERSE_HEXAGON is the rotated form, at 90 degrees).  Since the goal here is "packing", the choice should clearly be for HEXAGON (since they are the bestagons ) And how do you go from hexagon to circle? Ahh, well there are a couple of issues to consider: Regular hexagons are composed of six equilateral triangles, so "side" and "radius" are actually the same The radius of a hexagon is the radius of the circumscribing circle, not the inscribed circle, and circumscribed circles of neighboring tiles would overlap (doh!) The Generate Tessellation tool uses hexagon area to determine hexagon size, not side length/radius* *This is cool, though, because there's a formula to derive area ( A ) from side ( a ) So what's the key thing to know? In order to use hexagons to generate packed circles, we need to define the hexagon with respect to the midpoint of the hexagon side (where the inscribed circle touches).  And what is the difference? I'll save the math, for now, and cut to the chase: You need to increase the circle radius by the inverse of the sine of the interior angle size of a hexagon -- 60 degrees (aka pi / 3 radians).  In Python, that works out to:   import math diameter = 100 radius = diameter * 0.5 side = radius / math.sin(math.pi / 3.0) hex_area = 1.5 * math.sqrt(3) * side * side     The rest of the problem is simple enough to script: Generate Tessellation with the compensated area Extract centroids from the hexagons as points Use the points to buffer circles of the original diameter And that would look like this (with an arbitrary Cartesian CRS and extent, and kludged naming):   sr = arcpy.SpatialReference(3968) tess_shape = 'TRANSVERSE_HEXAGON' # or 'HEXAGON' arcpy.env.workspace = r"C:\Temp\packed_circles\cartesian.gdb" extent = '-1000 -1000 1000 1000' subscript = 1 print("Generating hexagons with side {:.6f}...".format(side)) arcpy.management.GenerateTessellation( Output_Feature_Class="hexagons{:d}".format(subscript), Extent=extent, Shape_Type=tess_shape, Size="{:f} SquareMeters".format(hex_area), Spatial_Reference=sr.exportToString() ) print("Generating centroids...") arcpy.management.FeatureToPoint( in_features="hexagons{:d}".format(subscript), out_feature_class="centroids{:d}".format(subscript), point_location="CENTROID" ) print("Generating circles (radius {:.1f})...".format(radius)) arcpy.analysis.Buffer( in_features="centroids{:d}".format(subscript), out_feature_class="circles{:d}".format(subscript), buffer_distance_or_field="{:.1f} Meters".format(radius), line_side="FULL", line_end_type="ROUND", dissolve_option="NONE", dissolve_field=None, method="PLANAR" )   Which gives us 100m diameter circles, packed tightly: And how do I use these circles to find the maximum number of circles in an arbitrary polygon?  Well, that's going to be a bit trickier (aka, "I haven't coded it yet"), so I'll address it in Part II. - V ... View more