Circle packing
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
class Point:
def __init__(self, x=0, y=0):
self.x = x
self.y = y
def squareDistTo(self, x, y):
return (self.x-x)**2 + (self.y-y)**2
class PackingCircle:
def __init__(self, x, y, r=1):
self.x = x
self.y = y
self.r = r
def __contains__(self, point):
return point.squareDistTo(self.x, self.y) < self.r**2
def plot_circle(x, y, r, fig):
circle=plt.Circle((x, y), r, fill=False)
fig.gca().add_artist(circle)
def plot_packingcircle(circle, fig):
circle=plt.Circle((circle.x, circle.y), circle.r, fill=False)
fig.gca().add_artist(circle)
def plot_regularpolygon(x, y, numVert, r, fig):
polygon= mpatches.RegularPolygon( (x, y), numVert, r, fill=True, alpha=0.2)
fig.gca().add_artist(polygon)
def plot_hexagonwithcircle(x, y, r, ol, fig):
plot_regularpolygon(x, y, 6, r*2*ol, fig)
a = PackingCircle (x, y, r)
b = PackingCircle (0+x, y+2*r * ol, r)
c = PackingCircle (0+x, y-2*r * ol, r)
d = PackingCircle (x+r * math.sqrt(3) *ol, y+r*ol, r)
e = PackingCircle (x+r * math.sqrt(3)*ol, y-r*ol, r)
f = PackingCircle (x-r * math.sqrt(3)*ol, y+r*ol, r)
g = PackingCircle (x-r * math.sqrt(3)*ol, y-r*ol, r)
plot_packingcircle(a, fig)
plot_packingcircle(b, fig)
plot_packingcircle(c, fig)
plot_packingcircle(d, fig)
plot_packingcircle(e, fig)
plot_packingcircle(f, fig)
plot_packingcircle(g, fig)
return [a, b, c, d, e, f, g], r*2*ol
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[-5,5],ylim=[-5,5])
# plot_circle(0,0, 1, fig)
# plot_regularpolygon(2, 2, 6, 3, fig)
circles, polygon_r = plot_hexagonwithcircle(2, 2, 1, 0.8, fig)
print polygon_r
1.6
# intersection area calculation
# verify the correctness
import math
def floatRange(start, stop, numofsteps):
i = float(start)
step = (float(stop)-start)/numofsteps
while i < stop:
yield i
i += step
def calculateArea(circles, xlim, ylim):
inSelf = inNone = 0.0
for i in floatRange(0, xlim, 1000):
for j in floatRange(0, ylim, 1000):
p = Point(i, j)
for packingcicle in circles:
# if p.squareDistTo(packingcicle.x, packingcicle.y) < packingcicle.r**2:
if p in packingcicle:
inSelf += 1
break
else:
inNone += 1
return xlim * ylim * inSelf / (inSelf+inNone)
x, y, u, w = [1, 1, 2, 2]
a = PackingCircle(x,y, 1)
b = PackingCircle(u,w, 1)
print calculateArea([a,b], 5, 5)
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[0,5],ylim=[0,5])
plot_circle(1, 1, 1, fig)
plot_circle(2, 2, 1, fig)
5.69962006026
# fig.gca(xlim=[0,100],ylim=[0,100])
# verify it with hexagonal packing
import math
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[-50,50],ylim=[-50,50])
a = PackingCircle (0, 0, 10)
b = PackingCircle (0, 20, 10)
c = PackingCircle (0, -20, 10)
d = PackingCircle (10 * math.sqrt(3), 10, 10)
e = PackingCircle (10 * math.sqrt(3), -10, 10)
f = PackingCircle (-10 * math.sqrt(3), 10, 10)
g = PackingCircle (-10 * math.sqrt(3), -10, 10)
# b = PackingCircle (20, 0, 10)
# c = PackingCircle (-20, 0, 10)
# d = PackingCircle (10, 10 * math.sqrt(3), 10)
# e = PackingCircle (-10, 10 * math.sqrt(3), 10)
# f = PackingCircle (10, -10 * math.sqrt(3), 10)
# g = PackingCircle (-10, -10 * math.sqrt(3), 10)
#
plot_packingcircle(a, fig)
plot_packingcircle(b, fig)
plot_packingcircle(c, fig)
plot_packingcircle(d, fig)
plot_packingcircle(e, fig)
plot_packingcircle(f, fig)
plot_packingcircle(g, fig)
plot_regularpolygon(0, 0, 6, 20, fig)
polygon= mpatches.RegularPolygon( (0, 0), 6, 20, fill=True, alpha=0.2)
# get a overlapped hexagon
# assume unit hexagon input
# reference http://stackoverflow.com/questions/5193331/is-a-point-inside-regular-hexagon
def ifInsideHexagon(x, y, radius = 1):
x = x/radius
y = y/radius
# Check length (squared) against inner and outer radius
l2 = x * x + y * y
if (l2 > 1.0): return False
if (l2 < 0.75): return True
#Check against borders
px = x * 1.15470053838;
if (px > 1.0 or px < -1.0): return False
py = 0.5 * px + y;
if (py > 1.0 or py < -1.0): return False
if (px - py > 1.0 or px - py < -1.0): return False
return True
#verify correctness
# print ifInsideHexagon(0, 1)
# print ifInsideHexagon(0, -1)
# print ifInsideHexagon(0, 1.00001)
print ifInsideHexagon(math.sqrt(3)/2*2, 0.5*2, 2)
print ifInsideHexagon(math.sqrt(3)/2-0.001, 0.49999)
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[-1.5,1.5],ylim=[-1.5,1.5])
# plot_packingcircle(g, fig)
plot_regularpolygon(0, 0, 6, 1, fig)
plot_circle(math.sqrt(3)/2, 0.45, 0.01, fig)
False True
#measure the coverage.
def measureCoverage(circles, hexagon_r, xmin, xmax, ymin, ymax):
inCircle = inPoly = 0.0
for i in floatRange(xmin, xmax, 1000):
for j in floatRange(ymin, ymax, 1000):
p = Point(i, j)
if ifInsideHexagon(p.x, p.y, hexagon_r):
for packingcicle in circles:
if p in packingcicle:
inCircle += 1
break
else:
inPoly += 1
return inCircle / (inCircle+inPoly)
r= 1
overlap = 0.1
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[-5,5],ylim=[-5,5])
#without overlap
a = PackingCircle (0, 0, r)
b = PackingCircle (0, 2*r, r)
c = PackingCircle (0, -2*r, r)
d = PackingCircle (r * math.sqrt(3), r, r)
e = PackingCircle (r * math.sqrt(3), -r, r)
f = PackingCircle (-r * math.sqrt(3), r, r)
g = PackingCircle (-r * math.sqrt(3), -r, r)
plot_packingcircle(a, fig)
plot_packingcircle(b, fig)
plot_packingcircle(c, fig)
plot_packingcircle(d, fig)
plot_packingcircle(e, fig)
plot_packingcircle(f, fig)
plot_packingcircle(g, fig)
plot_regularpolygon(0, 0, 6, r*2, fig)
print measureCoverage([a,b,c,d,e,f,g], r*2, -5, 5, -5, 5)
0.906835700967
r= 1
ol = 0.86
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[-5,5],ylim=[-5,5])
#with overlap
a = PackingCircle (0, 0, r)
b = PackingCircle (0, 2*r * ol, r)
c = PackingCircle (0, -2*r * ol, r)
d = PackingCircle (r * math.sqrt(3) *ol, r*ol, r)
e = PackingCircle (r * math.sqrt(3)*ol, -r*ol, r)
f = PackingCircle (-r * math.sqrt(3)*ol, r*ol, r)
g = PackingCircle (-r * math.sqrt(3)*ol, -r*ol, r)
plot_packingcircle(a, fig)
plot_packingcircle(b, fig)
plot_packingcircle(c, fig)
plot_packingcircle(d, fig)
plot_packingcircle(e, fig)
plot_packingcircle(f, fig)
plot_packingcircle(g, fig)
plot_regularpolygon(0, 0, 6, r*2*ol, fig)
print measureCoverage([a,b,c,d,e,f,g], r*2*ol, -5, 5, -5, 5)
1.0
r= 1
ol = 0.96
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[-5,5],ylim=[-5,5])
#with overlap
a = PackingCircle (0, 0, r)
b = PackingCircle (0, 2*r * ol, r)
c = PackingCircle (0, -2*r * ol, r)
d = PackingCircle (r * math.sqrt(3) *ol, r*ol, r)
e = PackingCircle (r * math.sqrt(3)*ol, -r*ol, r)
f = PackingCircle (-r * math.sqrt(3)*ol, r*ol, r)
g = PackingCircle (-r * math.sqrt(3)*ol, -r*ol, r)
plot_packingcircle(a, fig)
plot_packingcircle(b, fig)
plot_packingcircle(c, fig)
plot_packingcircle(d, fig)
plot_packingcircle(e, fig)
plot_packingcircle(f, fig)
plot_packingcircle(g, fig)
plot_regularpolygon(0, 0, 6, r*2*ol, fig)
print measureCoverage([a,b,c,d,e,f,g], r*2*ol, -5, 5, -5, 5)
0.955883580221
r= 1
ol = 0.96
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[-5,5],ylim=[-5,5])
circles, polygon_r = plot_hexagonwithcircle(0, 0, 1, ol, fig)
circles_1_total = []
circles_2_total = []
for circle in circles:
circles_1, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_1_total.extend(circles_1)
for circle in circles_1_total:
circles_2, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_2_total.extend(circles_2)
# plot_regularpolygon(0, 0, 6, r*2*ol, fig)
# print measureCoverage([a,b,c,d,e,f,g], r*2*ol, -5, 5, -5, 5)
r= 1
ol = 0.96
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[-5,5],ylim=[-5,5])
circles, polygon_r = plot_hexagonwithcircle(-4, -4, 1, ol, fig)
circles_1_total = []
circles_2_total = []
for circle in circles:
circles_1, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_1_total.extend(circles_1)
for circle in circles_1_total:
circles_2, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_2_total.extend(circles_2)
r= 1
ol = 0.96
fig = plt.gcf()
fig.set_size_inches(5, 5)
fig.gca(xlim=[-5,5],ylim=[-5,5])
circles, polygon_r = plot_hexagonwithcircle(-4-0.3, -4, 1, ol, fig)
circles_1_total = []
circles_2_total = []
circles_3_total = []
circles_4_total = []
circles_5_total = []
circles_6_total = []
circles_7_total = []
circles_8_total = []
for circle in circles:
if circle.x < -4 or circle.y < -4 or circle.x > 5 or circle.y >5:
continue
circles_1, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_1_total.extend(circles_1)
print len(circles_1_total)
for circle in circles_1_total:
if circle.x < -4 or circle.y < -4 or circle.x > 5 or circle.y >5:
continue
circles_2, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_2_total.extend(circles_2)
print len(circles_2_total)
for circle in circles_2_total:
if circle.x < -4 or circle.y < -4 or circle.x > 5 or circle.y >5:
continue
circles_3, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_3_total.extend(circles_3)
print len(circles_3_total)
for circle in circles_3_total:
if circle.x < -4 or circle.y < -4 or circle.x > 5 or circle.y >5:
continue
circles_4, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_4_total.extend(circles_4)
print len(circles_4_total)
for circle in circles_4_total:
if circle.x < -4 or circle.y < -4 or circle.x > 5 or circle.y >5:
continue
circles_5, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_5_total.extend(circles_5)
print len(circles_4_total)
for circle in circles_5_total:
if circle.x < -4 or circle.y < -4 or circle.x > 5 or circle.y >5:
continue
circles_6, polygon_r = plot_hexagonwithcircle(circle.x, circle.y, 1, ol, fig)
circles_6_total.extend(circles_6)
7 28 140 763
3*math.sqrt(3)/2*(1.92*10)**2
957.7548145532863
100*100/957.7548145532863
10.441085597323962
10.441085597323962 *3
31.323256791971886