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 *


    

    




Populating the interactive namespace from numpy and matplotlib

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]:





array([  3,   6,  14,  16,  18,  19,  20,  22,  32,  33,  34,  38,  41,
        43,  57,  58,  70,  73,  82,  83, 104, 106, 109, 112])

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]:





<matplotlib.collections.PathCollection at 0x9263400>






    

In [12]:

    

scatter(width, width/length)


    

    
Out[12]:





<matplotlib.collections.PathCollection at 0x960c780>






    

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]:





[<matplotlib.lines.Line2D at 0x951c1d0>]






    

In [14]:

    

minVal = 3
maxVal = 15
numZones = 1000


    

In [15]:

    

# 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[15]:





42.196132780708453

In [16]:

    

cutout[0].area


    

    
Out[16]:





40.994999999999976

In [17]:

    

plot(zoneWidth,zoneWeight)


    

    
Out[17]:





[<matplotlib.lines.Line2D at 0x9620ba8>]






    

In [18]:

    

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 [19]:

    

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 [20]:

    

middle_points = zeros([len(cutoutRef),2])

for j in range(len(cutoutRef)):
    
    i = cutoutRef[j]
    
    print("Cutout %.f" % (i))
    output = optimise_middle(cutout[i],3)
    middle_points[j,:] = output.x
    print(" ")


    

    




Cutout 3
first local minima of -30.3181 at:
[a,b] = [0.2030, -0.1662]
agreeing local minima of -30.3181 at:
[a,b] = [0.2030, -0.1662]
agreeing local minima of -30.3181 at:
[a,b] = [0.2030, -0.1662]
 
Cutout 6
first local minima of -45.6274 at:
[a,b] = [0.4048, -0.5731]
agreeing local minima of -45.6274 at:
[a,b] = [0.4049, -0.5732]
agreeing local minima of -45.6274 at:
[a,b] = [0.4049, -0.5732]
 
Cutout 14
first local minima of -27.4354 at:
[a,b] = [0.2154, -0.4741]
agreeing local minima of -27.4354 at:
[a,b] = [0.2154, -0.4741]
agreeing local minima of -27.4354 at:
[a,b] = [0.2154, -0.4741]
 
Cutout 16
first local minima of -43.7230 at:
[a,b] = [0.0550, -0.6293]
agreeing local minima of -43.7230 at:
[a,b] = [0.0550, -0.6293]
agreeing local minima of -43.7230 at:
[a,b] = [0.0550, -0.6293]
 
Cutout 18
first local minima of -46.3919 at:
[a,b] = [-0.0911, -0.4511]
agreeing local minima of -46.3919 at:
[a,b] = [-0.0911, -0.4511]
agreeing local minima of -46.3919 at:
[a,b] = [-0.0911, -0.4511]
 
Cutout 19
first local minima of -41.8733 at:
[a,b] = [0.5059, 0.2835]
agreeing local minima of -41.8733 at:
[a,b] = [0.5059, 0.2835]
agreeing local minima of -41.8733 at:
[a,b] = [0.5059, 0.2835]
 
Cutout 20
first local minima of -25.8190 at:
[a,b] = [-0.1952, -0.3417]
agreeing local minima of -25.8190 at:
[a,b] = [-0.1952, -0.3417]
agreeing local minima of -25.8190 at:
[a,b] = [-0.1952, -0.3417]
 
Cutout 22
first local minima of -39.7169 at:
[a,b] = [-0.8338, -0.3404]
agreeing local minima of -39.7169 at:
[a,b] = [-0.8338, -0.3404]
agreeing local minima of -39.7169 at:
[a,b] = [-0.8338, -0.3404]
 
Cutout 32
first local minima of -46.6405 at:
[a,b] = [0.7673, 0.2137]
agreeing local minima of -46.6405 at:
[a,b] = [0.7673, 0.2137]
agreeing local minima of -46.6405 at:
[a,b] = [0.7673, 0.2137]
 
Cutout 33
first local minima of -37.2148 at:
[a,b] = [-0.4707, 0.3438]
agreeing local minima of -37.2148 at:
[a,b] = [-0.4707, 0.3438]
agreeing local minima of -37.2148 at:
[a,b] = [-0.4707, 0.3438]
 
Cutout 34
first local minima of -45.2489 at:
[a,b] = [0.3041, -0.3796]
agreeing local minima of -45.2489 at:
[a,b] = [0.3041, -0.3796]
agreeing local minima of -45.2489 at:
[a,b] = [0.3041, -0.3796]
 
Cutout 38
first local minima of -29.9645 at:
[a,b] = [-0.0210, -0.1653]
agreeing local minima of -29.9645 at:
[a,b] = [-0.0210, -0.1653]
agreeing local minima of -29.9645 at:
[a,b] = [-0.0210, -0.1653]
 
Cutout 41
first local minima of -45.3349 at:
[a,b] = [0.0229, -0.0468]
agreeing local minima of -45.3349 at:
[a,b] = [0.0229, -0.0468]
agreeing local minima of -45.3349 at:
[a,b] = [0.0229, -0.0468]
 
Cutout 43
first local minima of -41.8051 at:
[a,b] = [-0.1068, 0.1081]
agreeing local minima of -41.8051 at:
[a,b] = [-0.1068, 0.1082]
agreeing local minima of -41.8051 at:
[a,b] = [-0.1068, 0.1081]
 
Cutout 57
first local minima of -20.7966 at:
[a,b] = [-0.2085, -0.1790]
agreeing local minima of -20.7966 at:
[a,b] = [-0.2085, -0.1790]
agreeing local minima of -20.7966 at:
[a,b] = [-0.2085, -0.1790]
 
Cutout 58
first local minima of -40.1427 at:
[a,b] = [-0.0756, 0.1855]
agreeing local minima of -40.1427 at:
[a,b] = [-0.0756, 0.1855]
agreeing local minima of -40.1427 at:
[a,b] = [-0.0756, 0.1855]
 
Cutout 70
first local minima of -46.0033 at:
[a,b] = [-0.0318, -0.0744]
agreeing local minima of -46.0033 at:
[a,b] = [-0.0318, -0.0744]
agreeing local minima of -46.0033 at:
[a,b] = [-0.0318, -0.0744]
 
Cutout 73
first local minima of -45.9770 at:
[a,b] = [-0.2967, -0.3890]
agreeing local minima of -45.9770 at:
[a,b] = [-0.2967, -0.3890]
agreeing local minima of -45.9770 at:
[a,b] = [-0.2967, -0.3890]
 
Cutout 82
first local minima of -42.6682 at:
[a,b] = [-0.4789, -0.4374]
agreeing local minima of -42.6682 at:
[a,b] = [-0.4789, -0.4374]
agreeing local minima of -42.6682 at:
[a,b] = [-0.4789, -0.4374]
 
Cutout 83
first local minima of -40.7417 at:
[a,b] = [-0.9655, -0.3627]
agreeing local minima of -40.7417 at:
[a,b] = [-0.9655, -0.3627]
agreeing local minima of -40.7417 at:
[a,b] = [-0.9655, -0.3627]
 
Cutout 104
first local minima of -40.4247 at:
[a,b] = [-0.7267, 0.5328]
agreeing local minima of -40.4247 at:
[a,b] = [-0.7267, 0.5328]
agreeing local minima of -40.4247 at:
[a,b] = [-0.7267, 0.5328]
 
Cutout 106
first local minima of -43.1247 at:
[a,b] = [-0.9173, -0.3189]
agreeing local minima of -43.1247 at:
[a,b] = [-0.9173, -0.3189]
agreeing local minima of -43.1247 at:
[a,b] = [-0.9172, -0.3189]
 
Cutout 109
first local minima of -42.6864 at:
[a,b] = [0.3963, -0.0756]
agreeing local minima of -42.6864 at:
[a,b] = [0.3963, -0.0756]
agreeing local minima of -42.6864 at:
[a,b] = [0.3963, -0.0756]
 
Cutout 112
first local minima of -18.6816 at:
[a,b] = [0.2003, -0.0606]
agreeing local minima of -18.6816 at:
[a,b] = [0.2003, -0.0606]
agreeing local minima of -18.6816 at:
[a,b] = [0.2003, -0.0606]
 

In [21]:

    

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],linewidth=2)

plot(a,b,'go')

scaled_fig_end(10,10)


    

In [22]:

    

i


    

    
Out[22]:





106

In [23]:

    

testcutout = af.translate(cutout[i], xoff=-a, yoff=-b)
testcutout


    

    
Out[23]:

In [24]:

    

ratioAreaInZones = zeros(numZones)

for k in range(numZones):
    
    ratioAreaInZones[k] = testcutout.intersection(zone[k]).area / zone[k].area

plot(zoneWidth,ratioAreaInZones)


    

    
Out[24]:





[<matplotlib.lines.Line2D at 0x96a43c8>]






    

In [25]:

    

dotAngle = 90-90*ratioAreaInZones
plot(zoneWidth,dotAngle)


    

    
Out[25]:





[<matplotlib.lines.Line2D at 0x96edef0>]






    

In [26]:

    

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[26]:





(-4.0, 4.0, -8.0, 8.0)






    

In [27]:

    

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 [28]:

    

straightCutout = shapely_cutout(ordered_x,ordered_y)
# straightCutout = af.rotate(straightCutout, -45)
straightCutout


    

    
Out[28]:

In [29]:

    

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[29]:





 minimization_failures: 0
                   nit: 3
                     x: array(-0.4466198956424787)
               message: ['requested number of basinhopping iterations completed successfully']
                  nfev: 93
                   fun: array(-42.458294041863)

In [30]:

    

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 [31]:

    

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)
scaled_fig_end(10,10)


    

In [32]:

    

centred_cutout = [0.]*len(cutoutRef)
straightened_cutout = [0.]*len(cutoutRef)

visual_angle = zeros(len(cutoutRef))

for j in range(len(cutoutRef)):
    
    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])


    

    




Cutout 3
Cutout 6
Cutout 14
Cutout 16
Cutout 18
Cutout 19
Cutout 20
Cutout 22
Cutout 32
Cutout 33
Cutout 34
Cutout 38
Cutout 41
Cutout 43
Cutout 57
Cutout 58
Cutout 70
Cutout 73
Cutout 82
Cutout 83
Cutout 104
Cutout 106
Cutout 109
Cutout 112

In [33]:

    

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 [34]:

    

# zoneWeight = 3000*widthFit.derivative()(zoneWidth)**2
# plot(zoneWidth,zoneWeight)


    

In [35]:

    

# 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])



# 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 [36]:

    

ellipse_results = [0,]*len(cutoutRef)

for j in range(len(cutoutRef)):
    i = cutoutRef[j]
    
    print("Cutout %.f" % (i))
    ellipse_results[j] = straight_ellipse_calc(straightened_cutout[j],1,2,visual_angle[j])
    print(" ")


    

    




Cutout 3
first local minima of 61.9533 at:
[a,b] = [0.0000, 0.0000], [width, length] = [4.4987, 6.0771], theta = 171.8465

agreeing local minima of 61.9533 at:
[a,b] = [0.0000, 0.0000], [width, length] = [4.4987, 6.0771], theta = 171.8465

 
Cutout 6
first local minima of 556.2569 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.9962, 10.6694], theta = 17.8357

agreeing local minima of 556.2569 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.9962, 10.6694], theta = 17.8357

 
Cutout 14
first local minima of 62.8098 at:
[a,b] = [0.0000, 0.0000], [width, length] = [4.3887, 5.2880], theta = 142.2589

agreeing local minima of 62.8098 at:
[a,b] = [0.0000, 0.0000], [width, length] = [4.3887, 5.2880], theta = 142.2589

 
Cutout 16
first local minima of 253.5485 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.1078, 10.5159], theta = 170.9191

agreeing local minima of 253.5485 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.1078, 10.5159], theta = 170.9191

 
Cutout 18
first local minima of 180.1784 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.6026, 10.0758], theta = 43.6695

agreeing local minima of 180.1784 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.6026, 10.0758], theta = 43.6695

 
Cutout 19
first local minima of 237.8095 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.7998, 10.3750], theta = 86.0336

agreeing local minima of 237.8095 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.7998, 10.3750], theta = 86.0336

 
Cutout 20
first local minima of 43.6950 at:
[a,b] = [0.0000, 0.0000], [width, length] = [4.2890, 5.1124], theta = 36.0110

agreeing local minima of 43.6950 at:
[a,b] = [0.0000, 0.0000], [width, length] = [4.2890, 5.1124], theta = 36.0110

 
Cutout 22
first local minima of 342.8536 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.1781, 10.2609], theta = 118.5026

agreeing local minima of 342.8536 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.1781, 10.2609], theta = 118.5026

 
Cutout 32
first local minima of 940.1030 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.7152, 11.1657], theta = 65.2545

agreeing local minima of 940.1030 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.7152, 11.1657], theta = 65.2545

 
Cutout 33
first local minima of 27.9050 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.9057, 6.2095], theta = 5.9845

agreeing local minima of 27.9050 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.9057, 6.2095], theta = 5.9845

 
Cutout 34
first local minima of 246.9361 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.0447, 9.6741], theta = 49.5904

agreeing local minima of 246.9361 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.0447, 9.6741], theta = 49.5904

 
Cutout 38
first local minima of 61.3259 at:
[a,b] = [0.0000, 0.0000], [width, length] = [4.6775, 5.7242], theta = 115.0504

agreeing local minima of 61.3259 at:
[a,b] = [0.0000, 0.0000], [width, length] = [4.6775, 5.7242], theta = 115.0504

 
Cutout 41
first local minima of 101.4130 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.9191, 10.0118], theta = 146.9746

agreeing local minima of 101.4130 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.9191, 10.0118], theta = 146.9746

 
Cutout 43
first local minima of 43.0571 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.1002, 7.7357], theta = 115.6672

agreeing local minima of 43.0571 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.1002, 7.7357], theta = 115.6672

 
Cutout 57
first local minima of 42.8046 at:
[a,b] = [0.0000, 0.0000], [width, length] = [3.5747, 4.6188], theta = 0.6015

agreeing local minima of 42.8046 at:
[a,b] = [0.0000, 0.0000], [width, length] = [3.5747, 4.6188], theta = 0.6015

 
Cutout 58
first local minima of 95.4860 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.8804, 7.4549], theta = 29.8452

agreeing local minima of 95.4860 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.8804, 7.4549], theta = 29.8452

 
Cutout 70
first local minima of 68.1437 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.4610, 8.3632], theta = 169.8589

agreeing local minima of 68.1437 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.4610, 8.3632], theta = 169.8589

 
Cutout 73
first local minima of 265.4707 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.2941, 9.1378], theta = 109.0199

agreeing local minima of 265.4707 at:
[a,b] = [0.0000, 0.0000], [width, length] = [7.2941, 9.1378], theta = 109.0199

 
Cutout 82
first local minima of 130.8533 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.2688, 7.9747], theta = 177.1163

agreeing local minima of 130.8533 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.2688, 7.9747], theta = 177.1163

 
Cutout 83
first local minima of 110.5096 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.7823, 7.9311], theta = 165.7175

agreeing local minima of 110.5096 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.7823, 7.9311], theta = 165.7175

 
Cutout 104
first local minima of 153.2798 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.6589, 8.3283], theta = 155.0258

agreeing local minima of 153.2798 at:
[a,b] = [0.0000, 0.0000], [width, length] = [5.6589, 8.3283], theta = 155.0258

 
Cutout 106
first local minima of 435.2211 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.0872, 9.7082], theta = 154.4106

agreeing local minima of 435.2211 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.0872, 9.7082], theta = 154.4106

 
Cutout 109
first local minima of 127.2634 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.3068, 7.8272], theta = 38.3637

agreeing local minima of 127.2634 at:
[a,b] = [0.0000, 0.0000], [width, length] = [6.3068, 7.8272], theta = 38.3637

 
Cutout 112
first local minima of 16.2029 at:
[a,b] = [0.0000, 0.0000], [width, length] = [3.5387, 4.1496], theta = 100.8248

agreeing local minima of 16.2029 at:
[a,b] = [0.0000, 0.0000], [width, length] = [3.5387, 4.1496], theta = 100.8248

 

In [37]:

    

ellipse_def = zeros([len(cutoutRef),5])

for j in range(len(cutoutRef)):
    ellipse_def[j,:] = array(ellipse_results[j].x)


    

In [38]:

    

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],'r--')
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1],'g')

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 [39]:

    

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 [40]:

    

width[0:24] - new_width


    

    
Out[40]:





array([-0.01432214, -0.147601  ,  0.00259141,  0.14693614, -0.13047532,
       -0.25075809, -0.02924766, -0.06806439, -0.18854346, -0.00087269,
       -0.13177121, -0.02467374, -0.10171873, -0.00491839, -0.00985364,
        0.00935551, -0.08438723, -0.1036374 ,  0.0445798 , -0.11523122,
       -0.21834959, -0.17281856,  0.01468404, -0.00077132])

In [41]:

    

cutoutRef


    

    
Out[41]:





array([  3,   6,  14,  16,  18,  19,  20,  22,  32,  33,  34,  38,  41,
        43,  57,  58,  70,  73,  82,  83, 104, 106, 109, 112])

In [42]:

    

# 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


    

    
Out[42]:






  

    

      

      
cutoutRef
      
width
      
new_width
      
length
      
new_length
      
factor
    
  
  

    

      
0 
      
   3
      
 4.484404
      
 4.498726
      
  6.136510
      
  6.077114
      
 1.032704
    
    

      
1 
      
   6
      
 6.848634
      
 6.996235
      
 10.332762
      
 10.669437
      
 1.001302
    
    

      
2 
      
  14
      
 4.391292
      
 4.388701
      
  5.270661
      
  5.288022
      
 1.043027
    
    

      
3 
      
  16
      
 6.254735
      
 6.107799
      
 10.368288
      
 10.515871
      
 0.999567
    
    

      
4 
      
  18
      
 7.472170
      
 7.602645
      
 10.011654
      
 10.075798
      
 0.998054
    
    

      
5 
      
  19
      
 5.549091
      
 5.799849
      
 10.134021
      
 10.375006
      
 1.001089
    
    

      
6 
      
  20
      
 4.259799
      
 4.289047
      
  5.122864
      
  5.112375
      
 1.037967
    
    

      
7 
      
  22
      
 5.110011
      
 5.178076
      
 10.251878
      
 10.260920
      
 1.011713
    
    

      
8 
      
  32
      
 7.526681
      
 7.715225
      
 10.998181
      
 11.165712
      
 0.993330
    
    

      
9 
      
  33
      
 5.904805
      
 5.905678
      
  6.208632
      
  6.209543
      
 1.006789
    
    

      
10
      
  34
      
 6.912947
      
 7.044718
      
  9.578994
      
  9.674088
      
 0.993330
    
    

      
11
      
  38
      
 4.652777
      
 4.677451
      
  5.727794
      
  5.724189
      
 1.027190
    
    

      
12
      
  41
      
 6.817334
      
 6.919052
      
 10.014400
      
 10.011766
      
 0.994829
    
    

      
13
      
  43
      
 6.095283
      
 6.100202
      
  7.754348
      
  7.735733
      
 1.007783
    
    

      
14
      
  57
      
 3.564813
      
 3.574667
      
  4.626695
      
  4.618824
      
 1.059704
    
    

      
15
      
  58
      
 5.889774
      
 5.880419
      
  7.441537
      
  7.454923
      
 1.008999
    
    

      
16
      
  70
      
 7.376619
      
 7.461006
      
  8.625289
      
  8.363232
      
 0.997190
    
    

      
17
      
  73
      
 7.190506
      
 7.294143
      
  9.602273
      
  9.137754
      
 0.995472
    
    

      
18
      
  82
      
 6.313361
      
 6.268781
      
  7.884968
      
  7.974677
      
 1.000867
    
    

      
19
      
  83
      
 5.667032
      
 5.782263
      
  8.415000
      
  7.931109
      
 1.010552
    
    

      
20
      
 104
      
 5.440583
      
 5.658933
      
  9.067747
      
  8.328343
      
 1.006601
    
    

      
21
      
 106
      
 5.914386
      
 6.087204
      
  9.902046
      
  9.708158
      
 1.000871
    
    

      
22
      
 109
      
 6.321438
      
 6.306754
      
  7.786481
      
  7.827168
      
 1.006979
    
    

      
23
      
 112
      
 3.537958
      
 3.538729
      
  4.146843
      
  4.149551
      
 1.053904
    
  

In [42]: