In [1]:
%pylab inline
import pandas as pd
import csv
import shapely.affinity as af
import shapely.geometry as sh
from scipy.interpolate import UnivariateSpline
from scipy.optimize import basinhopping
# from equivalent_ellipse import *
from scaled_figures import *
In [35]:
pylab.rcParams['savefig.dpi'] = 254
In [2]:
def shapely_cutout(x,y):
cutoutPolyList = []
for i in range(len(x)):
cutoutPolyList = cutoutPolyList + [(x[i],y[i])]
cutout = sh.Polygon(cutoutPolyList)
return cutout
unit_circle = sh.Point(0,0).buffer(1)
def create_ellipse(x):
a = x[0]
b = x[1]
width = x[2]
length = x[3]
theta = x[4]
ellipse = af.scale(unit_circle, xfact = width/2, yfact = length/2)
ellipse = af.translate(ellipse, xoff=a, yoff=b)
ellipse = af.rotate(ellipse, theta)
return ellipse
In [3]:
x_list = list()
y_list = list()
with open('../data/custom_cutout_x.csv', 'r') as x_csvfile:
with open('../data/custom_cutout_y.csv', 'r') as y_csvfile:
x_reader = csv.reader(x_csvfile, delimiter=',', lineterminator='\n')
y_reader = csv.reader(y_csvfile, delimiter=',', lineterminator='\n')
for row in x_reader:
x_list += [row]
for row in y_reader:
y_list += [row]
# num_cutouts = len(x_list)
num_cutouts = 1
x_array = [0,]*num_cutouts
y_array = [0,]*num_cutouts
for i in range(num_cutouts):
x_array[i] = array(x_list[i], dtype='float')
y_array[i] = array(y_list[i], dtype='float')
cutout = [0,]*num_cutouts
for i in range(num_cutouts):
cutout[i] = shapely_cutout(x_array[i],y_array[i])
In [4]:
# Use pandas to load the csv file
cutoutData = pd.read_csv("cutoutdata.csv")
# Define width, length, and factor as the csv columns
# which had that corresponding header
width = cutoutData['width'].values
length = cutoutData['length'].values
factor = cutoutData['factor'].values
# Swap width and length values if order was mistaken
swap = width > length
store = width[swap]
width[swap] = length[swap]
length[swap] = store
In [5]:
scatter(width,factor)
widthFit = UnivariateSpline(width,factor)
width_i = linspace(2,10)
widthFit_i = widthFit(width_i)
plot(width_i,widthFit_i,'r-')
Out[5]:
In [6]:
test = widthFit.derivative()
plot(width_i,test(width_i)**2)
Out[6]:
In [7]:
# Create Zones
minVal = 1
maxVal = sqrt(2) * 10 * 2
numZones = 200
zoneWidth = linspace(minVal,maxVal,numZones)
zoneWeight = widthFit.derivative()(zoneWidth)**2
circle = [0,]*(numZones)
circle[0] = sh.Point(0,0).buffer(minVal/2)
zone = [0,]*numZones
zone[0] = circle[0]
for i in range(1,numZones-1):
circle[i] = sh.Point(0,0).buffer(zoneWidth[i]/2)
zone[i] = circle[i].difference(circle[i-1])
zone[-1] = sh.Point(0,0).buffer(20).difference(circle[-2])
weightedArea = 0.
for i in range(numZones):
weightedArea += cutout[0].intersection(zone[i]).area * zoneWeight[i]
weightedArea
Out[7]:
In [8]:
zone[28]
Out[8]:
In [9]:
def weightedAreaCalc(inputShape):
weightedArea = 0.
for i in range(numZones):
weightedArea += inputShape.intersection(zone[i]).area * zoneWeight[i]
return weightedArea
def weightedAreaDifference(cutout,ellipse,a,b):
cutout = af.translate(cutout, xoff=-a, yoff=-b)
ellipse = af.translate(ellipse, xoff=-a, yoff=-b)
return (weightedAreaCalc(ellipse.difference(cutout)) +
weightedAreaCalc(cutout.difference(ellipse)))
In [10]:
cutout[0].area
Out[10]:
In [10]:
In [11]:
def ellipse_in_cutout(minimiserInput,cutout,weight):
ellipse = create_ellipse(minimiserInput)
a = minimiserInput[0]
b = minimiserInput[1]
toMinimise = 100*weightedAreaDifference(cutout,ellipse,a,b)
return toMinimise
class MyTakeStep(object):
def __init__(self, stepsize=0.5):
self.stepsize = stepsize
def __call__(self, x):
s = self.stepsize
x[0] += np.random.uniform(-3*s, 3*s)
x[1] += np.random.uniform(-3*s, 3*s)
x[2] += np.random.uniform(-6*s, 6*s)
x[3] += np.random.uniform(-8*s, 8*s)
x[4] += np.random.uniform(-90*s, 90*s)
return x
def print_fun(x, f, accepted):
global successCount_basinhopping
global minVals_basinhopping
global nSuccess_basinhopping
a = x[0]
b = x[1]
w = abs(x[2])
l = abs(x[3])
if (l > w):
width = w
length = l
theta = np.mod(x[4],180)
else:
width = l
length = w
theta = np.mod(x[4]+90,180)
f = np.around(f,4)
if accepted:
if np.all(np.isnan(minVals_basinhopping)):
minVals_basinhopping[0] = f
successCount_basinhopping = 1
print(("first local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
elif (np.nanmin(minVals_basinhopping) == f):
minVals_basinhopping[successCount_basinhopping] = f
successCount_basinhopping += 1
print(("agreeing local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
elif (np.nanmin(minVals_basinhopping) > f):
minVals_basinhopping[0] = f
successCount_basinhopping = 1
print(("new local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
else:
print(("rejected local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
else:
print(("failed to find local minima at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (a,b,width,length,theta))
if successCount_basinhopping >= nSuccess_basinhopping:
return True
def eq_ellipse_calc(cutout,weight,userNSuccess):
mytakestep = MyTakeStep()
global nSuccess_basinhopping
nSuccess_basinhopping = userNSuccess
global minVals_basinhopping
minVals_basinhopping = np.empty(nSuccess_basinhopping)
minVals_basinhopping[:] = np.nan
global successCount_basinhopping
successCount_basinhopping = 0
minimizer_kwargs = {"args": (cutout,weight), "method": 'L-BFGS-B', "bounds": ((-5,5),(-5,5),(2,10),(2.5,10),(None,None)) }
x0 = np.array([0,0,3,5,0])
output = basinhopping(ellipse_in_cutout,x0,
niter=1000,minimizer_kwargs=minimizer_kwargs,
take_step=mytakestep, callback=print_fun)
return output
In [12]:
ellipse_results = [0,]*num_cutouts
for i in range(num_cutouts):
print("Cutout %.f" % (i))
ellipse_results[i] = eq_ellipse_calc(cutout[i],1,5)
print(" ")
In [13]:
ellipse_def = zeros([num_cutouts,5])
for i in range(num_cutouts):
ellipse_def[i,:] = array(ellipse_results[i].x)
In [36]:
i = 0
ellipse = create_ellipse(ellipse_def[i,:])
scaled_fig_start(10,10)
plot(cutout[i].exterior.xy[0],cutout[i].exterior.xy[1])
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1],'r--')
text(-4.5,4,('Width: %.2f cm' % (ellipse_def[i,2])),fontsize=11)
text(-4.5,4.5,('Length: %.2f cm' % (ellipse_def[i,3])),fontsize=11)
plot(ellipse.centroid.xy[0],ellipse.centroid.xy[1],'go')
scaled_fig_end(10,10)
savefig('../figures/straightened/initial.png')
In [15]:
a = ellipse_def[:,0]
b = ellipse_def[:,1]
width = zeros([num_cutouts])
length = zeros([num_cutouts])
theta = zeros([num_cutouts])
w = abs(ellipse_def[:,2])
l = abs(ellipse_def[:,3])
ref = l > w
width[ref] = w[ref]
length[ref] = l[ref]
theta[ref] = mod(ellipse_def[:,4],180)[ref]
width[~ref] = l[~ref]
length[~ref] = w[~ref]
theta[~ref] = mod(ellipse_def[:,4]+90,180)[~ref]
In [16]:
testcutout = af.translate(cutout[i], xoff=-a[i], yoff=-b[i])
testellipse = af.translate(ellipse, xoff=-a[i], yoff=-b[i])
In [17]:
testcutout
Out[17]:
In [18]:
testcutout.intersection(zone[9])
Out[18]:
In [19]:
ratioAreaInZones = zeros(numZones)
for i in range(numZones):
ratioAreaInZones[i] = testcutout.intersection(zone[i]).area / zone[i].area
In [20]:
plot(zoneWidth,ratioAreaInZones)
Out[20]:
In [21]:
dotAngle = 90-90*ratioAreaInZones
plot(zoneWidth,dotAngle)
Out[21]:
In [22]:
x_store = array([])
y_store = array([])
for i in range(len(dotAngle)):
if (dotAngle[i] > 0) & (dotAngle[i] < 90):
hypot = zoneWidth[i]/2
y_val = hypot * sin(pi/180*dotAngle[i])
x_val = hypot * cos(pi/180*dotAngle[i])
x_store = append(x_store,x_val)
y_store = append(y_store,y_val)
x_mirror = array([x_val,x_val,-x_val,-x_val])
y_mirror = array([y_val,-y_val,y_val,-y_val])
scatter(x_mirror,y_mirror)
axis('equal')
Out[22]:
In [23]:
ordered_x = concatenate((x_store,-x_store[::-1],-x_store,x_store[::-1]))
ordered_y = concatenate((y_store,y_store[::-1],-y_store,-y_store[::-1]))
In [24]:
plot(ordered_x,ordered_y)
axis('equal')
Out[24]:
In [25]:
straightCutout = shapely_cutout(ordered_x,ordered_y)
straightCutout = af.rotate(straightCutout, -45)
straightCutout
Out[25]:
In [34]:
ratioAreaInZones_test = zeros(numZones)
for i in range(numZones):
ratioAreaInZones_test[i] = straightCutout.intersection(zone[i]).area / zone[i].area
plot(zoneWidth,ratioAreaInZones)
plot(zoneWidth,ratioAreaInZones_test,'x')
Out[34]:
In [26]:
def straight_ellipse_calc(cutout,weight,userNSuccess):
mytakestep = MyTakeStep()
global nSuccess_basinhopping
nSuccess_basinhopping = userNSuccess
global minVals_basinhopping
minVals_basinhopping = np.empty(nSuccess_basinhopping)
minVals_basinhopping[:] = np.nan
global successCount_basinhopping
successCount_basinhopping = 0
minimizer_kwargs = {"args": (cutout,weight), "method": 'L-BFGS-B', "bounds": ((0,0),(0,0),(2,15),(2.5,15),(None,None)) }
x0 = np.array([0,0,3,5,0])
output = basinhopping(ellipse_in_cutout,x0,
niter=1000,minimizer_kwargs=minimizer_kwargs,
take_step=mytakestep, callback=print_fun)
return output
In [27]:
output = straight_ellipse_calc(straightCutout,1,5)
output
Out[27]:
In [37]:
ellipse = create_ellipse(output.x)
scaled_fig_start(10,10)
plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1])
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1],'r--')
text(-4,3,('Width: %.2f cm' % (output.x[2])),fontsize=11)
text(-4,3.5,('Length: %.2f cm' % (output.x[3])),fontsize=11)
plot(ellipse.centroid.xy[0],ellipse.centroid.xy[1],'go')
scaled_fig_end(10,10)
savefig('../figures/straightened/final.png')
In [29]:
# a = ellipse_def[:,0]
# b = ellipse_def[:,1]
# width = zeros([num_cutouts])
# length = zeros([num_cutouts])
# theta = zeros([num_cutouts])
# w = abs(ellipse_def[:,2])
# l = abs(ellipse_def[:,3])
# ref = l > w
# width[ref] = w[ref]
# length[ref] = l[ref]
# theta[ref] = mod(ellipse_def[:,4],180)[ref]
# width[~ref] = l[~ref]
# length[~ref] = w[~ref]
# theta[~ref] = mod(ellipse_def[:,4]+90,180)[~ref]
In [30]:
# box_bounds = zeros([num_cutouts])
# for i in range(num_cutouts):
# box_bounds[i] = max(
# array(
# [cutout[i].bounds[2]-cutout[i].bounds[0],
# cutout[i].bounds[3]-cutout[i].bounds[1]]
# )
# )
In [31]:
# cutout_dimensions = pd.DataFrame(
# array([a,b,width,length,theta,box_bounds]).transpose(),
# columns=['a','b','width','length','theta','box_bounds']
# )
In [32]:
# cutout_dimensions.to_csv('../data/custom_cutout_dimensions.csv')
In [32]: