Circle packing

In [134]:
%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
In [98]:
# 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
In [100]:
# 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)
In [108]:
# 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

stacking with overlap

In [110]:
#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)
In [111]:
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
In [124]:
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
In [118]:
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

stacking with boundary

In [141]:
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)
In [142]:
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)
In [156]:
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
In [159]:
3*math.sqrt(3)/2*(1.92*10)**2
Out[159]:
957.7548145532863
In [160]:
100*100/957.7548145532863
Out[160]:
10.441085597323962
In [161]:
10.441085597323962 *3
Out[161]:
31.323256791971886