In [1]:
%pylab inline
from scipy.optimize import basinhopping
from scipy.interpolate import UnivariateSpline
import pandas as pd
import csv
import shapely.affinity as af
import shapely.geometry as sh
from equivalent_ellipse import *
from scaled_figures import *
from descartes.patch import PolygonPatch
In [2]:
pylab.rcParams['savefig.dpi'] = 110
In [3]:
x_list = list()
y_list = list()
with open('../data/cutout_x.csv', 'r') as x_csvfile:
with open('../data/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)
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
cutoutRef = cutoutData['Reference'].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]:
cutoutRef = cutoutRef[~isnan(cutoutRef)].astype(int)
cutoutRef
Out[5]:
In [6]:
# factor
In [7]:
# with open('cutout_factors','r') as f:
# loaded_factor = eval(f.read())
# dictArray = array(list(loaded_factor.items()),float)
# cutoutRef = dictArray[:,0].astype(int)
# factor = dictArray[:,1]
# # width = width[cutoutRef]
# # length = length[cutoutRef]
# factor
In [8]:
# with open('custom_cutout_factors','r') as f:
# loaded_custom_factor = eval(f.read())
# customDimensionsDict = {'5x13':[5,13],
# '5x10':[5,10],
# '5x8':[5,8],
# '4x13':[4,13],
# '4x10':[4,10],
# '4x8':[4,8],
# '4x6.5':[4,6.5],
# '3x13':[3,13],
# '3x9':[3,9],
# '3x6.5':[3,6.5],
# '3x5':[3,5]}
# # customDimensionsDict = {'5x8':[5,8],
# # '4x8':[4,8],
# # '4x6.5':[4,6.5],
# # '3x9':[3,9],
# # '3x6.5':[3,6.5],
# # '3x5':[3,5]}
# customDataLabelList = list(customDimensionsDict.keys())
# custom_factors = [loaded_custom_factor[x] for x in customDataLabelList]
# custom_widths = [customDimensionsDict[x][0] for x in customDataLabelList]
# custom_lengths = [customDimensionsDict[x][1] for x in customDataLabelList]
# ellipse_cutouts = [0,]*len(custom_widths)
# for i in range(len(custom_widths)):
# ellipse_cutouts[i] = create_ellipse([0,0,custom_widths[i],custom_lengths[i],-45])
# cutoutRef = append(cutoutRef,arange(len(cutout),len(cutout)+len(ellipse_cutouts)))
# cutout.extend(ellipse_cutouts)
# factor = append(factor,custom_factors)
# width = append(width,custom_widths)
# length = append(length,custom_lengths)
In [9]:
# testCutout = cutout[70]
# testCutout
In [10]:
with open('circle_cutout_factors','r') as f:
loaded_circle_factor = eval(f.read())
circleDictArray = array(list(loaded_circle_factor.items()),float)
circleDiameter = circleDictArray[:,0]
circleFactor = circleDictArray[:,1]
circle_cutouts = [0,]*len(circleDiameter)
for i in range(len(circleDiameter)):
circle_cutouts[i] = create_ellipse([0,0,circleDiameter[i],circleDiameter[i],0])
In [11]:
scatter(width,length)
Out[11]:
In [12]:
scatter(width, width/length)
Out[12]:
In [13]:
scatter(width, 1/factor)
circleFit = UnivariateSpline(circleDiameter,1/circleFactor)
diameter_i = linspace(circleDiameter.min(),circleDiameter.max())
plot(diameter_i,circleFit(diameter_i),'r--')
widthFit = UnivariateSpline(width,1/factor)
width_i = linspace(width.min(),width.max())
plot(width_i,widthFit(width_i),'k--')
Out[13]:
In [14]:
pylab.rcParams['savefig.dpi'] = 90
In [15]:
ratio = width/length
artificial = isnan(cutoutData['Reference'].values)
scatter(width[~artificial],ratio[~artificial])
circleRef = ratio == 1
scatter(width[circleRef],ratio[circleRef],c='red')
scatter(width[artificial & ~circleRef],ratio[artificial & ~circleRef],c='green')
title("Distrbution of selected cutouts")
xlabel("Width (cm)")
ylabel("Ratio")
legend(['Clinical','Circles','Ellipses'],loc=4)
Out[15]:
In [16]:
scatter(width[~artificial],1/factor[~artificial])
circleRef = ratio == 1
scatter(width[circleRef],1/factor[circleRef],c='red')
scatter(width[artificial & ~circleRef],1/factor[artificial & ~circleRef],c='green')
# widthFit = UnivariateSpline(width,1/factor)
# width_i = linspace(width.min(),width.max())
plot(width_i,widthFit(width_i),'r--')
title("The effect of width on cutout factor")
ylabel("1 / cutout factor")
xlabel("Width cm")
Out[16]:
In [77]:
# scatter(width, 1/factor)
pylab.rcParams['savefig.dpi'] = 90
# widthFit = UnivariateSpline(width,1/factor)
width_i = linspace(width.min(),width.max(),1000)
derivFunc = widthFit.derivative()
deriv = widthFit.derivative()(width_i)
deriv[deriv<0] = 0
plot(width_i,deriv,'b')
plot(width_i,zeros(len(width_i)),'k--')
# scatter(width_i,deriv,c=deriv,edgecolors='none',s=2)
title("Contribution to cutout factor\nfrom ring of given diameter")
ylabel("Derivative of (1 / cutout factor wrt width)")
xlabel("Diameter | Width (cm)")
xlim([3,9])
ylim([-0.002,0.025])
Out[77]:
In [17]:
In [18]:
minVal = 3
maxVal = 15
numZones = 1000
In [19]:
# Create Zones
zoneWidth = linspace(minVal,maxVal,numZones)
zoneWeight = 70*widthFit.derivative()(zoneWidth)
zoneWeight[zoneWeight < 0] = 0
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[19]:
In [20]:
cutout[0].area
Out[20]:
In [21]:
plot(zoneWidth,zoneWeight)
Out[21]:
In [22]:
def weightedAreaOptimise(x,inputShape):
a = x[0]
b = x[1]
movedShape =af.translate(inputShape, xoff=-a, yoff=-b)
weightedArea = 0.
for i in range(numZones):
weightedArea += movedShape.intersection(zone[i]).area * zoneWeight[i]
return -weightedArea
In [23]:
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)
return x
def print_fun(x, f, accepted):
global successCount_basinhopping
global minVals_basinhopping
global nSuccess_basinhopping
a = x[0]
b = x[1]
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]")
% (f,a,b))
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]")
% (f,a,b))
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]")
% (f,a,b))
else:
print(("rejected local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f]")
% (f,a,b))
else:
print(("failed to find local minima at:\n"+
"[a,b] = [%.4f, %.4f]")
% (a,b))
if successCount_basinhopping >= nSuccess_basinhopping:
return True
def optimise_middle(cutout,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,), "method": 'L-BFGS-B',
"bounds": ((-5,5),(-5,5)) }
x0 = np.array([0,0])
output = basinhopping(weightedAreaOptimise,x0,
niter=1000,minimizer_kwargs=minimizer_kwargs,
take_step=mytakestep, callback=print_fun)
return output
In [24]:
middle_points = zeros([len(cutoutRef),2])
# for j in range(len(cutoutRef)):
for j in [21]:
i = cutoutRef[j]
print("Cutout %.f" % (i))
output = optimise_middle(cutout[i],3)
middle_points[j,:] = output.x
print(" ")
In [54]:
j = 21
i = cutoutRef[j]
a,b = middle_points[j,:]
scaled_fig_start(10,10)
plot(cutout[i].exterior.xy[0],cutout[i].exterior.xy[1])
plot(a,b,'go')
scaled_fig_end(10,10)
In [26]:
testcutout = af.translate(cutout[i], xoff=-a, yoff=-b)
testcutout
Out[26]:
In [27]:
ratioAreaInZones = zeros(numZones)
for k in range(numZones):
ratioAreaInZones[k] = testcutout.intersection(zone[k]).area / zone[k].area
plot(zoneWidth,ratioAreaInZones)
Out[27]:
In [28]:
dotAngle = 90-90*ratioAreaInZones
plot(zoneWidth,dotAngle)
Out[28]:
In [29]:
x_store = array([])
y_store = array([])
for k in range(len(dotAngle)):
if (dotAngle[k] > 0) & (dotAngle[k] < 90):
hypot = zoneWidth[k]/2
y_val = hypot * sin(pi/180*dotAngle[k])
x_val = hypot * cos(pi/180*dotAngle[k])
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[29]:
In [30]:
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 [31]:
straightCutout = shapely_cutout(ordered_x,ordered_y)
# straightCutout = af.rotate(straightCutout, -45)
straightCutout
Out[31]:
In [32]:
def visualRotation(x,original,new):
rotatedShape = af.rotate(new, x)
shapeIntersection = rotatedShape.intersection(original)
# weightedOverlap = 0.
# for i in range(numZones):
# weightedOverlap += shapeIntersection.intersection(zone[i]).area * zoneWeight[i]
# print(-shapeIntersection.area)
return -shapeIntersection.area
minimizer_kwargs = {"args": (testcutout,straightCutout,), "method": 'Powell'}
x0 = np.array([-90])
output = basinhopping(visualRotation,x0,
niter=3,minimizer_kwargs=minimizer_kwargs)
straightCutout = af.rotate(straightCutout, output.x)
output
Out[32]:
In [33]:
scaled_fig_start(10,10)
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],linewidth=2)
plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=1)
scatter(0,0,c='green',s=30)
scaled_fig_end(10,10)
In [34]:
scaled_fig_start(10,10)
# plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=0.5)
scatter(0,0,c='green',s=40)
for k in range(0,numZones,100):
plot(zone[k].exterior.xy[0],zone[k].exterior.xy[1],'k',linewidth=2,alpha=0.5)
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],linewidth=2)
overlap = zone[2]
scaled_fig_end(10,10)
In [35]:
# scatter(width, 1/factor)
# widthFit = UnivariateSpline(width,1/factor)
width_i = linspace(width.min(),15,1000)
derivFunc = widthFit.derivative()
deriv = widthFit.derivative()(width_i)
# plot(width_i,deriv,'b')
plot(width_i,zeros(len(width_i)),'k--')
scatter(width_i,deriv,c=deriv,edgecolors='none',s=2)
title("Contribution to cutout factor\nfrom ring of given diameter")
ylabel("Derivative of (1 / cutout factor wrt width)")
xlabel("Diameter | Width (cm)")
# xlim([3,15])
# ylim([-0.005,0.025])
Out[35]:
In [36]:
numRands = 100000
x0 = 0
y0 = 0
dist = width.min()/2 + rand(numRands) * width.ptp()/2
angle = 2*pi*rand(numRands)
xVals = x0 + dist*cos(angle)
yVals = y0 + dist*sin(angle)
scaled_fig_start(10,10)
scatter(xVals,yVals,c=derivFunc(dist),edgecolors='none',s=5,alpha=0.1)
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],'k',linewidth=2)
scaled_fig_end(10,10)
In [36]:
In [38]:
scaled_fig_start(10,10)
# plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=0.5)
scatter(0,0,c='green',s=40)
for k in range(0,numZones,100):
plot(zone[k].exterior.xy[0],zone[k].exterior.xy[1],'k',linewidth=1,alpha=0.5)
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],linewidth=2)
overlap = zone[2]
scaled_fig_end(10,10)
In [39]:
numZonesTest = 10
zoneWidthTest = linspace(minVal,maxVal,numZonesTest)
# zoneWeight = 70*widthFit.derivative()(zoneWidth)
circleTest = [0,]*(10)
circleTest[0] = sh.Point(0,0).buffer(minVal/2)
zoneTest = [0,]*numZonesTest
zoneTest[0] = circleTest[0]
for i in range(1,numZonesTest-1):
circleTest[i] = sh.Point(0,0).buffer(zoneWidthTest[i]/2)
zoneTest[i] = circleTest[i].difference(circleTest[i-1])
zoneTest[-1] = sh.Point(0,0).buffer(20).difference(circleTest[-2])
In [40]:
pylab.rcParams['savefig.dpi'] = 70
BLUE = '#6699cc'
YELLOW = '#ffcc33'
GREEN = '#339933'
GRAY = '#999999'
RED = '#ff3333'
In [41]:
scaled_fig_start(10,10)
ax = gca()
# plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=0.5)
scatter(0,0,c='green',s=40)
for k in range(0,numZonesTest):
plot(zoneTest[k].exterior.xy[0],zoneTest[k].exterior.xy[1],'k',linewidth=1,alpha=0.5)
intersect = testcutout.intersection(zoneTest[2])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=GREEN, ec=GREEN, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=GREEN, ec=GREEN, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
intersect = testcutout.intersection(zoneTest[3])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
intersect = testcutout.intersection(zoneTest[5])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=RED, ec=RED, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=RED, ec=RED, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],linewidth=2)
scaled_fig_end(10,10)
In [42]:
scaled_fig_start(10,10)
ax = gca()
# plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=0.5)
scatter(0,0,c='green',s=40)
for k in range(0,numZonesTest):
plot(zoneTest[k].exterior.xy[0],zoneTest[k].exterior.xy[1],'k',linewidth=1,alpha=0.5)
intersect = straightCutout.intersection(zoneTest[2])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=GREEN, ec=GREEN, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=GREEN, ec=GREEN, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
intersect = straightCutout.intersection(zoneTest[3])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
intersect = straightCutout.intersection(zoneTest[5])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=RED, ec=RED, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=RED, ec=RED, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],linewidth=2)
scaled_fig_end(10,10)
In [43]:
centred_cutout = [0.]*len(cutoutRef)
straightened_cutout = [0.]*len(cutoutRef)
visual_angle = zeros(len(cutoutRef))
# for j in range(len(cutoutRef)):
for j in [21]:
i = cutoutRef[j]
print("Cutout %d" %i)
a,b = middle_points[j,:]
centred_cutout[j] = af.translate(cutout[i], xoff=-a, yoff=-b)
ratioAreaInZones = zeros(numZones)
for k in range(numZones):
ratioAreaInZones[k] = centred_cutout[j].intersection(zone[k]).area / zone[k].area
dotAngle = 90-90*ratioAreaInZones
x_store = array([])
y_store = array([])
for k in range(len(dotAngle)):
if (dotAngle[k] > 0) & (dotAngle[k] < 90):
hypot = zoneWidth[k]/2
y_val = hypot * sin(pi/180*dotAngle[k])
x_val = hypot * cos(pi/180*dotAngle[k])
x_store = append(x_store,x_val)
y_store = append(y_store,y_val)
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]))
straightened_cutout[j] = shapely_cutout(ordered_x,ordered_y)
minimizer_kwargs = {"args": (centred_cutout[j],straightened_cutout[j],), "method": 'Powell'}
x0 = np.array([-90])
output = basinhopping(visualRotation,x0,
niter=3,minimizer_kwargs=minimizer_kwargs)
visual_angle[j] = output.x
straightened_cutout[j] = af.rotate(straightened_cutout[j], visual_angle[j])
In [44]:
# j = find(cutoutRef == 38)[0]
# i = cutoutRef[j]
# plot(centred_cutout[j].exterior.xy[0],centred_cutout[j].exterior.xy[1])
# plot(straightened_cutout[j].exterior.xy[0],straightened_cutout[j].exterior.xy[1],'r--')
# axis('equal')
# grid(True)
In [45]:
# zoneWeight = 3000*widthFit.derivative()(zoneWidth)**2
# plot(zoneWidth,zoneWeight)
In [60]:
# numZones = 1000
# zoneWidth = linspace(minVal,maxVal,numZones)
# zoneWeight = zoneWeight = 3000*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])
zoneWidth = linspace(minVal,maxVal,numZones)
zoneWeight = 70*widthFit.derivative()(zoneWidth)
zoneWeight[zoneWeight < 0] = 0
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])
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)))
def weightedAreaDifference(cutout,ellipse,a,b):
cutout = af.translate(cutout, xoff=-a, yoff=-b)
ellipse = af.translate(ellipse, xoff=-a, yoff=-b)
return ellipse.difference(cutout).area + cutout.difference(ellipse).area
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 straight_ellipse_calc(cutout,weight,userNSuccess,angle):
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),
(minVal,maxVal),
(minVal,maxVal),
(angle,angle)) }
x0 = np.array([0,0,3,5,angle])
output = basinhopping(ellipse_in_cutout,x0,
niter=1000,minimizer_kwargs=minimizer_kwargs,
take_step=mytakestep, callback=print_fun)
return output
In [61]:
ellipse_results = [0,]*len(cutoutRef)
# for j in range(len(cutoutRef)):
for j in [21]:
i = cutoutRef[j]
print("Cutout %.f" % (i))
ellipse_results[j] = straight_ellipse_calc(straightened_cutout[j],1,2,visual_angle[j])
print(" ")
In [62]:
ellipse_def = zeros([len(cutoutRef),5])
# for j in range(len(cutoutRef)):
for j in [21]:
ellipse_def[j,:] = array(ellipse_results[j].x)
In [71]:
# j = find(cutoutRef == 14)[0]
# i = cutoutRef[j]
ellipse = create_ellipse(ellipse_def[j,:])
scaled_fig_start(12,12)
# plot(centred_cutout[j].exterior.xy[0],centred_cutout[j].exterior.xy[1])
plot(straightened_cutout[j].exterior.xy[0],straightened_cutout[j].exterior.xy[1],'b')
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1],'g',linewidth=4,alpha=0.5)
text(-5,5,('Width: %.2f cm' % (ellipse_def[j,2])),fontsize=11)
text(-5,4.5,('Length: %.2f cm' % (ellipse_def[j,3])),fontsize=11)
scaled_fig_end(12,12)
In [53]:
# a = ellipse_def[:,0]
# b = ellipse_def[:,1]
# new_width = zeros([len(cutoutRef)])
# new_length = zeros([len(cutoutRef)])
# theta = zeros([len(cutoutRef)])
# w = abs(ellipse_def[:,2])
# l = abs(ellipse_def[:,3])
# ref = l > w
# new_width[ref] = w[ref]
# new_length[ref] = l[ref]
# theta[ref] = mod(ellipse_def[:,4],180)[ref]
# new_width[~ref] = l[~ref]
# new_length[~ref] = w[~ref]
# theta[~ref] = mod(ellipse_def[:,4]+90,180)[~ref]
In [ ]:
# width[0:24] - new_width
In [ ]:
# cutoutRef
In [ ]:
# # Display table
# df = pd.DataFrame(transpose([cutoutRef, width[0:24], new_width, length[0:24], new_length, factor[0:24]]))
# df.columns = ['cutoutRef','width', 'new_width', 'length', 'new_length' , 'factor']
# df.to_csv(path_or_buf = "new_cutoutdata.csv", index=False)
# df
In [ ]: