In [1]:
%pylab inline
import csv, os
from glob import glob
import pandas as pd
import shapely.affinity as af
import shapely.geometry as sh
from equivalent_ellipse import *
from scaled_figures import *
# pylab.rcParams['savefig.dpi'] = 254
pylab.rcParams['savefig.dpi'] = 100
In [2]:
from descartes.patch import PolygonPatch
from scipy.spatial import Delaunay
from scipy.interpolate import InterpolatedUnivariateSpline
In [3]:
dimensions = array([6,10,14,20,25])
cornerRadii = array([0.5,0.5,0.5,0.5,0.5])
cornerType = "round" # Options: round, mitred, none
In [4]:
pos = dimensions[0]/2 - cornerRadii[0]
example = sh.Polygon([(-pos, -pos),
(pos, -pos),
(pos, pos),
(-pos, pos)])
example = example.buffer(cornerRadii[0], join_style=3)
example
Out[4]:
In [5]:
scaled_fig_start(7,7)
plot(example.exterior.xy[0],example.exterior.xy[1])
scaled_fig_end(7,7)
In [5]:
In [5]:
In [5]:
In [6]:
def ellipse_in_applicator(minimiserInput,applicator,weight):
ellipse = create_ellipse(minimiserInput)
return weight*ellipse.difference(applicator).area + applicator.difference(ellipse).area
class MyTakeStep(object):
def __init__(self, stepsize=0.5):
self.stepsize = stepsize
def __call__(self, x):
s = self.stepsize
x[0] += np.random.uniform(-0*s, 0*s)
x[1] += np.random.uniform(-0*s, 0*s)
x[2] += np.random.uniform(-0*s, 0*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 optimise_largest_insert(applicator,weight,userNSuccess,width):
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": (applicator,weight), "method": 'L-BFGS-B',
"bounds": ((0,0),(0,0),
(width,width),
(3,None),
(None, None)) }
x0 = np.array([0,0,width,5,0])
output = basinhopping(ellipse_in_applicator,x0,
niter=1000,minimizer_kwargs=minimizer_kwargs,
take_step=mytakestep, callback=print_fun)
return output
In [7]:
output = optimise_largest_insert(example,100,3,2)
output
Out[7]:
In [8]:
ellipse = create_ellipse(output.x)
ellipse
Out[8]:
In [9]:
scaled_fig_start(7,7)
plot(example.exterior.xy[0],example.exterior.xy[1])
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1])
scaled_fig_end(7,7)
In [10]:
output = optimise_largest_insert(example,100,3,5)
output
ellipse = create_ellipse(output.x)
ellipse
Out[10]:
In [11]:
scaled_fig_start(7,7)
plot(example.exterior.xy[0],example.exterior.xy[1])
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1])
scaled_fig_end(7,7)
In [11]:
In [12]:
erroded = example.buffer(-0.2)
erroded
Out[12]:
In [13]:
output = optimise_largest_insert(erroded,100,3,2)
output
Out[13]:
In [14]:
ellipse = create_ellipse(output.x)
ellipse
Out[14]:
In [15]:
scaled_fig_start(7,7)
plot(example.exterior.xy[0],example.exterior.xy[1])
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1])
plot(erroded.exterior.xy[0],erroded.exterior.xy[1])
scaled_fig_end(7,7)
In [16]:
maxThinRatio = output.x[2] / output.x[3]
maxThinLogRatio = log2(maxThinRatio)
maxThinLogRatio
Out[16]:
In [16]:
In [16]:
In [16]:
In [17]:
def ellipse_in_applicator_changeWidth(minimiserInput,applicator,weight):
adjInput = minimiserInput.copy()
adjInput[3] = minimiserInput[2]/minimiserInput[3]
ellipse = create_ellipse(adjInput)
return weight*ellipse.difference(applicator).area + applicator.difference(ellipse).area
class takeStep_changeWidth(object):
def __init__(self, stepsize=0.5):
self.stepsize = stepsize
def __call__(self, x):
s = self.stepsize
x[0] += np.random.uniform(-0*s, 0*s)
x[1] += np.random.uniform(-0*s, 0*s)
x[2] += np.random.uniform(-2*s, 2*s)
x[3] += np.random.uniform(-0*s, 0*s)
x[4] += np.random.uniform(-90*s, 90*s)
return x
def optimise_largest_insert_changeWidth(applicator,weight,userNSuccess,ratio):
mytakestep = takeStep_changeWidth()
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": (applicator,weight), "method": 'L-BFGS-B',
"bounds": ((0,0),(0,0),
(0.01,None),
(ratio,ratio),
(None, None)) }
x0 = np.array([0,0,3,ratio,0])
output = basinhopping(ellipse_in_applicator_changeWidth,x0,
niter=1000,minimizer_kwargs=minimizer_kwargs,
take_step=mytakestep, callback=print_fun)
return output
In [18]:
# output = optimise_largest_insert_changeWidth(paddedZone,100,3,2**mean(thinElipseLogRatios[-2::]))
# adjOutput = output.x.copy()
# adjOutput[3] = output.x[2]/output.x[3]
# ellipse = create_ellipse(adjOutput)
# ellipse
In [19]:
# output.x
In [20]:
# scaled_fig_start(appDimension+1,appDimension+1)
# plot(appature.exterior.xy[0],appature.exterior.xy[1])
# plot(paddedZone.exterior.xy[0],paddedZone.exterior.xy[1])
# plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1])
# scaled_fig_end(appDimension+1,appDimension+1)
In [21]:
def batch_max_ellipse(maximalEllpiseWidths,appature,appDimension):
maximalEllipseLengths = zeros(len(maximalEllpiseWidths))
for i in range(len(maximalEllpiseWidths)):
print("Maximal ellipse with width %f0" % (maximalEllpiseWidths[i]))
output = optimise_largest_insert(appature,100,3,maximalEllpiseWidths[i])
print(" ")
ellipse = create_ellipse(output.x)
scaled_fig_start(appDimension+1,appDimension+1)
plot(appature.exterior.xy[0],appature.exterior.xy[1])
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1])
scaled_fig_end(appDimension+1,appDimension+1)
show()
maximalEllipseLengths[i] = output.x[3]
return maximalEllipseLengths
In [22]:
def create_apature(appDimension):
bevel = 1.0
pos = appDimension/2 - bevel
appature = sh.Polygon([(-pos, -pos),
(pos, -pos),
(pos, pos),
(-pos, pos)])
appature = appature.buffer(bevel, join_style=3)
return appature
In [23]:
appDimension = 6
maximalEllpiseWidths6 = array([3,4,5])
appature = create_apature(appDimension)
maximalEllipseLengths6 = batch_max_ellipse(maximalEllpiseWidths6,
appature,
appDimension)
In [24]:
figure(figsize=(5, 5))
scatter(maximalEllpiseWidths6,maximalEllipseLengths6,s=70)
circles = array([3,4,5,6])
scatter(circles,circles,c='red',s=70)
# scatter(4,4,c='red')
scatter(3,5,c='green',s=70)
# scatter(15,21,c='green')
title('Initial measurement 6x6cm')
axis("equal")
xlabel("Width (cm)")
ylabel("Length (cm)")
print(maximalEllpiseWidths6)
print(maximalEllipseLengths6)
In [25]:
sqrt(2)*10 / 4
Out[25]:
In [26]:
appDimension = 10
maximalEllpiseWidths10 = array([4,6,8])
appature = create_apature(appDimension)
maximalEllipseLengths10 = batch_max_ellipse(maximalEllpiseWidths10,
appature,
appDimension)
In [27]:
figure(figsize=(5, 5))
# scatter(maximalEllpiseWidths6[[0,2]],maximalEllipseLengths6[[0,2]],c='yellow',s=70)
scatter([3,4.5],[7,6],c='yellow',s=70)
scatter(maximalEllpiseWidths10,maximalEllipseLengths10,s=70)
circles = array([6,7.5,9,10])
scatter(circles,circles,c='red',s=70)
scatter([3,3],[9,11],c='green',s=70)
title('Initial measurement 10x10cm')
axis("equal")
xlabel("Width (cm)")
ylabel("Length (cm)")
print(maximalEllpiseWidths10)
print(maximalEllipseLengths10)
In [28]:
sqrt(2)*15 / 4
Out[28]:
In [29]:
appDimension = 15
maximalEllpiseWidths15 = array([6,9,12])
appature = create_apature(appDimension)
maximalEllipseLengths15 = batch_max_ellipse(maximalEllpiseWidths15,
appature,
appDimension)
In [30]:
18/4
Out[30]:
In [31]:
figure(figsize=(5, 5))
# scatter(maximalEllpiseWidths10[[0,2]],maximalEllipseLengths10[[0,2]],c='yellow',s=70)
scatter([4,6.5],[12,10],c='yellow',s=70)
scatter(maximalEllpiseWidths15,maximalEllipseLengths15,c='blue',s=70)
circles = array([10,12,14,15])
scatter(circles,circles,c='red',s=70)
scatter([4.5,4],[18,15],c='green',s=70)
title('Initial measurement 15x15cm')
axis("equal")
xlabel("Width (cm)")
ylabel("Length (cm)")
print(maximalEllpiseWidths15)
print(maximalEllipseLengths15)
In [32]:
sqrt(2)*20 / 4
Out[32]:
In [33]:
appDimension = 20
maximalEllpiseWidths20 = array([7,12,16])
appature = create_apature(appDimension)
maximalEllipseLengths20 = batch_max_ellipse(maximalEllpiseWidths20,
appature,
appDimension)
In [34]:
6*4
Out[34]:
In [35]:
figure(figsize=(5, 5))
# scatter(maximalEllpiseWidths15[[0,2]],maximalEllipseLengths15[[0,2]],c='yellow',s=70)
scatter([6,10],[19,15],c='yellow',s=70)
scatter(maximalEllpiseWidths20,maximalEllipseLengths20, c='blue', s=70)
circles = array([15,18,20])
scatter(circles,circles,c='red',s=70)
scatter([6],[23],c='green',s=70)
title('Initial measurement 20x20cm')
axis("equal")
xlabel("Width (cm)")
ylabel("Length (cm)")
print(maximalEllpiseWidths20)
print(maximalEllipseLengths20)
In [36]:
sqrt(2)*25 / 4
Out[36]:
In [37]:
appDimension = 25
maximalEllpiseWidths25 = array([9,15,20])
appature = create_apature(appDimension)
maximalEllipseLengths25 = batch_max_ellipse(maximalEllpiseWidths25,
appature,
appDimension)
In [38]:
7.5 * 4
Out[38]:
In [39]:
figure(figsize=(5, 5))
# scatter(maximalEllpiseWidths20,maximalEllipseLengths20, c='yellow',s=70)
scatter([7,11,15],[27,23,20], c='yellow',s=70)
scatter(maximalEllpiseWidths25,maximalEllipseLengths25, c='blue', s=70)
circles = array([20,23,25])
scatter(circles,circles,c='red',s=70)
scatter([8],[31],c='green',s=70)
title('Initial measurement 25x25cm')
axis("equal")
xlabel("Width (cm)")
ylabel("Length (cm)")
print(maximalEllpiseWidths25)
print(maximalEllipseLengths25)
In [39]:
In [39]:
In [39]:
In [40]:
figure(figsize=(5, 5))
plot([3,3],[0,40],'k--')
plot([0,10],[0,40],'g')
plot([0,30],[0,30],'r')
# scatter(3,3,s=70)
scatter(maximalEllpiseWidths6,maximalEllipseLengths6,s=70)
scatter(maximalEllpiseWidths10,maximalEllipseLengths10,s=70)
scatter(maximalEllpiseWidths15,maximalEllipseLengths15,s=70)
scatter(maximalEllpiseWidths20,maximalEllipseLengths20,s=70)
scatter(maximalEllpiseWidths25,maximalEllipseLengths25,s=70)
circles = array([3,4,5,6,7.5,9,10,12,14,15,18,20,23,25])
# scatter()
# circles = array([20,23,25])
scatter(circles,circles,c='red',s=70)
# scatter([8],[31],c='green',s=70)
# title('Initial measurement 20x20cm')
scatter([3,3,3,4,4.5,6,8],[5,9,11,15,18,23,31],c='green',s=70)
# scatter([3,4.5,4,6.5,6,10,7,11,15],[7,6,12,10,19,15,27,23,20],c='yellow',s=70)
axis("equal")
xlabel("Width (cm)")
ylabel("Length (cm)")
xlim([0,30])
ylim([0,38])
title("Overview of desired clinical cutout regions")
Out[40]: