diff --git a/batch_coord_gen.py b/batch_coord_gen.py new file mode 100644 index 0000000..09760ed --- /dev/null +++ b/batch_coord_gen.py @@ -0,0 +1,289 @@ +''' +*batch_coord_gen.py* +Used to create an instance of a borefield + +Functions: + + genCases([float,float],[float,float],[float,float],[float,float],float,float, + float,float,float,float,float,float,float,float + ,float,float,float,float,float,float,string,string,[float,float] + ,[float,float],float,float,float,float,float,float, + float,float,float,float,float,string,boolean, + boolean,int,boolean,float,float,float) -> none + main() -> none + +''' +from math import pi +import math +import time +import coord_gen as InnerProgram +import csv +import random +from nogo import Shapes +from matplotlib.backends.backend_pdf import PdfPages + +def genCases (BaseSLocStart,BaseSLocStop,TopSLocStart,TopSLocStop,BaseWidthStart,BaseWidthStop, + TopWidthStart,TopWidthStop,YSpacStart,YSpacStop,XSpacStart,XSpacStop,xwStart,xwStop + ,ywStart,ywStop,cxStart,cxStop,cyStart,cyStop,BaseFileName="",Directory="",BaseSLocStep=[1,1] + ,TopSLocStep=[1,1],BaseWidthStep=1,TopWidthStep=1,YSpacStep=1,XSpacStep=1,xwStep=1,ywStep=1, + cxStep=1,cyStep=1,thetaStart=pi/8,thetaStop=pi/2.0,thetaStep=pi/11,unit="m",graph = False, + outputCases = True,SampleRate = 1000,Show = False,minArea = .45,maxArea=.55,mV = .75): + ''' + Function generates a series of x,y points repersenting a field of bore holes + in a trapezoidal shape. Returs empty if boreHole field does not meet given requirements + + Parameters + ------------- + :param BaseSLocStart: [float,float] + the initial value of the location of the base of the field + :param BaseSLocStop: [float,float] + the final value of the location of the base of the field + :param TopSLocStart: [float,float] + the initial value of the location of the top of the field + :param TopSLocStop: [float,float] + the final value of the location of the top of the field + :param BaseWidthStart: float + the initial value of the width of the base of the field + :param BaseWidthStop: float + the final value of the width of the base of the field + :param TopWidthStart: float + the initial value of the width of the top of the field + :param TopWidthStop: float + the final value fo the width of the top of the field + :param YSpacStart: float + the initial value for the minimum spacing in the y-dir + :param YSpacStop: float + the final value for the minimum spacing in the y-dir + :param XSpacStart: float + the intial value for the minimum spacing in the x-dir + :param XSpacStop: float + the final value for the minimum spacing int the x-dir + :param xwStart: float + the intial value for the xwidth of the nogo zone + :param xwStop: float + the final value for the xwidth of the nogo zone + :param ywStart: float + the initial value for the ywidth of the nogo zone + :param ywStop: float + the final value for the ywidth of the nogo zone + :param cxStart: float + the initial x value for the centroid of the nogo zone + :param cxStop: float + the final x value for the centroid of the nogo zone + :param cyStart: float + the initial y value for the centroid of the nogo zone + :param cyStop: float + the final y value for the centroid of the nogo zone + :param BaseFileName: string + the base filename to store bore fields (case number is appended to this) + :param Directory: string + the directory to place bore field csv files + :param BaseSLocStep: [float,float] + the stepsize to increment the x,y location of the base of the field + :param TopSLocStep: [float,float] + the stepsize to increment the x,y location of the top of the field + :param BaseWidthStep: float + the stepsize to increment the width of the base of the field + :param TopWidthStep: float + the stepsize to increment the width of the top of the field + :param YSpacStep: float + the stepsize to increment the y spacing of the field + :param XSpacStep: float + the stepsize to increment the x spacing of the field + :param xwStep: float + the stepsize to increment the xwidth of the nogo zone + :param ywStep: float + the stepsize to increment the ywidth of the nogo zone + :param cxStep: float + the stepsize to increment the x value of the centroid of the nogo zone + :param cyStep: float + the stepsize to increment the y vlaue of the centroid of the nogo zone + :param thetaStart: float + the initial value of the rotation of the nogo zone + :param thetaStop: float + the final value of the rotation of the nogo zone + :param thetaStep: float + the stepsize to increment the rotation of the nogo zone + :param unit: string + string containting the deisred unit only m is currently supported + :param graph: boolean + decides whether random cases should be selected and graphed for user inspection + :param outputCases: boolean + decides whether csv files should be written + :param SampleRate: int + integer controlling how often the generated cases are sampled + :param Show: boolean + decides whether graphed fields should be displayed (true) or thrown to file (false) + :param minArea: float + the minimum area of the field that the nogo zone should take up + :param maxArea: float + the maximum are of the field that the nogo zone should take up + :param mV: float + the fraction of vertices of the nogozone that should be inside the field + :return: nothing + ''' + + start = time.time() + BaseSLocStartx,BaseSLocStarty=BaseSLocStart + BaseSLocStopx,BaseSLocStopy = BaseSLocStop + TopSLocStartx,TopSLocStarty = TopSLocStart + TopSLocStopx,TopSLocStopy = TopSLocStop + BaseSLocStepx,BaseSLocStepy=BaseSLocStep + TopSLocStepx,TopSLocStepy = TopSLocStep + numberOfPos = ((math.floor((TopWidthStop-TopWidthStart)/TopWidthStep)+1)* (math.floor((BaseWidthStop-BaseWidthStart)/BaseWidthStep)+1) + * (math.floor((BaseSLocStopx-BaseSLocStartx)/BaseSLocStepx)+1) * (math.floor((BaseSLocStopy-BaseSLocStarty)/BaseSLocStepy)+1)* + (math.floor((TopSLocStopx-TopSLocStartx)/TopSLocStepx)+1)* (math.floor((TopSLocStopy-TopSLocStarty)/TopSLocStepy)+1)* + (math.floor((thetaStop - thetaStart) / thetaStep) + 1)*(math.floor((cxStop-cxStart)/cxStep)+1)*(math.floor((cyStop-cyStart)/cyStep)+1)* + (math.floor((ywStop - ywStart) / ywStep) + 1)*(math.floor((xwStop-xwStart)/xwStep)+1)*(math.floor((YSpacStop-YSpacStart)/YSpacStep)+1)* + (math.floor((XSpacStop - XSpacStart) / XSpacStep) + 1)*6 + ) + print("Number of Possible Fields",numberOfPos) + + + if unit == "m": + unit = "(m)" + else: + unit = "(ft)" + n = 0 + shapeStrings = ["S","U","T","L","Baseball","B"] + thrownOut = 0 + BoreHoleFields = {} + selected = [0] + selectCurrent = 0 + shapes = {} + thrownOut = 0 + smallArea=0 + + cy = cyStart + while cy <= cyStop: + cx = cxStart + while cx <= cxStop: + yw = ywStart + while yw <= ywStop: + xw = xwStart + while xw<=xwStop: + ang = thetaStart + while ang <= thetaStop: + for shapeString in shapeStrings: + shape = Shapes(cx, cy, xw, yw, theta=ang,sh=shapeString) + XSpac = XSpacStart + while XSpac <= XSpacStop: + YSpac = YSpacStart + while YSpac <= YSpacStop: + TopWidth = TopWidthStart + while TopWidth <= TopWidthStop: + BaseWidth = BaseWidthStart + while BaseWidth <= BaseWidthStop: + BaseSLocx = BaseSLocStartx + while BaseSLocx <= BaseSLocStopx: + BaseSLocy = BaseSLocStarty + while BaseSLocy <= BaseSLocStopy: + TopSLocx = TopSLocStartx + while TopSLocx <= TopSLocStopx: + TopSLocy = TopSLocStarty + while TopSLocy <= TopSLocStopy: + ar = shape.getArea()/( .5 * (BaseWidth + TopWidth) * ( + TopSLocy - BaseSLocy)) + if ar < minArea or ar > maxArea: + TopSLocy += TopSLocStepy + smallArea += 1 + continue + holes = InnerProgram.genBoreHoleConfig([BaseSLocx,BaseSLocy],BaseWidth,[TopSLocx,TopSLocy],TopWidth,YSpac,XSpac,shape,minVert=mV) #creates borefield + if len(holes)==0: + TopSLocy += TopSLocStepy + thrownOut+=1 + continue + if graph: + if n == selected[selectCurrent]: + shapes[selectCurrent] = shape + selected.append( + random.randrange(selected[selectCurrent]+1, selected[selectCurrent]+SampleRate, 1)) + selectCurrent+=1 + BoreHoleFields[n]=holes + n+=1 + TopSLocy+=TopSLocStepy + TopSLocx+=TopSLocStepx + BaseSLocy+=BaseSLocStepy + BaseSLocx+=BaseSLocStepx + BaseWidth+=BaseWidthStep + TopWidth+=TopWidthStep + YSpac += YSpacStep + XSpac+=XSpacStep + ang += thetaStep + xw+=xwStep + yw+=ywStep + cx += cxStep + cy += cyStep + print("Number of Cases thrown out: ",thrownOut) + print("Number with too small of area: ",smallArea) + print("Total Fields Generated: ", n) + print("Beginning Outputting") + if outputCases: + for i in range(n): + with open("".join([Directory, BaseFileName, "_", str(i), ".csv"]), "w") as outputFile: + csvWriter = csv.writer(outputFile) + csvWriter.writerow(["x" + unit, "y" + unit]) + csvWriter.writerows(BoreHoleFields[i]) + if graph: + selected.__delitem__(len(selected)-1) + if Show: + i=0 + for sel in selected: + print(sel) + InnerProgram.plotField(BoreHoleFields[sel], shapes[i], title="Case_" + str(sel)) + i+=1 + else: + with PdfPages(Directory + "Graphs" + ".pdf") as PdP: + i=0 + for sel in selected: + InnerProgram.plotField(BoreHoleFields[sel], shapes[i], SaveOrShow=Show, Pdf=PdP, title="Case_" + str(sel)) + i+=1 + end = time.time() + print("Time elapsed: ",end-start) + return +def main (): + ''' + An example use of the genCases function + :return: nothing + ''' + BSLSta = [0,0] + BSLSto = [0,0] + TSLSta= [0,60] + TSLSto= [40,100] + BWSta = 60 + BWSto = 100 + TWSta = 0 + TWSto = 100 + YSSta = 5 + YSSto = 6 + XSSta = YSSta + XSSto = YSSto + xwSta = 14 + xwSto = 30 + ywSta = 34 + ywSto = 50 + cxSta = 20 + cxSto = 80 + cySta = 20 + cySto = 80 + xwS = 10 + ywS = xwS + cxS = 30 + cyS = cxS + TWS = 50 + BWS = 40 + TopSLS = [40,40] + XSS = .5 + YSS = .5 + BFN = "\Files\Case_TTNGR" + BD = "D:\Work\GSHP\Optimizing\Cases" + genCases(BSLSta,BSLSto,TSLSta,TSLSto,BWSta,BWSto,TWSta,TWSto,YSSta,YSSto,XSSta,XSSto + ,xwSta,xwSto,ywSta,ywSto,cxSta,cxSto,cySta,cySto,graph = True, + BaseFileName = BFN,Directory = BD,xwStep = xwS,ywStep = ywS + ,cxStep = cxS, cyStep = cyS,outputCases=False,TopSLocStep = TopSLS + ,BaseWidthStep = BWS,TopWidthStep = TWS,SampleRate = 500,Show=True + ,XSpacStep=XSS,YSpacStep=YSS,minArea=.45,maxArea=.55,mV = .75) + + return +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/coord_gen.py b/coord_gen.py new file mode 100644 index 0000000..d8d2eed --- /dev/null +++ b/coord_gen.py @@ -0,0 +1,263 @@ +''' +*coord_gen.py* +Used to create an instance of a borefield + +Functions: + + genBoreHoleConfig(float,float,float,float,float,float,object,int,float) -> [[float]] + plotField([[float]],object,string,boolean,object) -> none + processRows(float,float,object,float,float,float,float,int) + def Between(float,float,float) -> double + distribute(float,float,float,float,[[float]]) -> appends to given array + +''' +import math +import numpy as np +from math import pi +import matplotlib.pyplot as plt +import operator +from nogo import Shapes +from matplotlib.backends.backend_pdf import PdfPages + +def genBoreHoleConfig(BaseSLoc,BaseWidth,TopSLoc,TopWidth,YSpac,XSpac,shape,maxLength = 250,minVert=.75): + ''' + Function generates a series of x,y points repersenting a field of bore holes + in a trapezoidal shape. Returs empty if boreHole field does not meet given requirements + + Parameters + ------------- + :param BaseSLoc: [float,float] + x,y coordinates of the leftmost point of the base of the trapezoid + :param BaseWidth: float + the width of the base of the trapezoid + :param TopSLoc: [float,float] + x,y coordinates of the leftmost point of the top of the trapezoid + :param TopWidth: float + the width of the top of the trapezoid + :param YSpac: float + the minimum spacing between points in the y-dir + :param XSpac: float + the minimum spacing between points in the x-dir + :param shape: shape object + the obstruction or "nogo zone" to the borefield + :param maxLength: int + the maximum number of boreholes allowed in the field + :param minVert: float + the fraction of vertices of the no-go zone required to be in the bore hole field + :return: [[float]] -> 2 col + n rows + ''' + minVert =minVert*len(shape.c) + + nrows = int((TopSLoc[1]-BaseSLoc[1])//YSpac) + rowspace = float((TopSLoc[1]-BaseSLoc[1]))/nrows + + s1 = ((TopSLoc[0]-BaseSLoc[0])/(TopSLoc[1]-BaseSLoc[1])) + s2 = ((TopSLoc[0]-BaseSLoc[0]+TopWidth-BaseWidth)/(TopSLoc[1]-BaseSLoc[1])) + nin = 0 + for cen in shape.c: + if cen[1] > TopSLoc[1] or cen[1] < BaseSLoc[1]: + continue + elif cen[0] > BaseSLoc[0] + s1*cen[1] and cen[0] < BaseSLoc[0]+BaseWidth + s2*cen[1]: + nin+=1 + if nin < minVert: + return [] + boreHoles =np.array( [row for i in range(nrows+1) if not (rows := ProcessRows(BaseSLoc,rowspace,shape,XSpac,s1,BaseWidth,s2,i))==None for row in rows]) + if len(boreHoles) > maxLength: + return [] + return boreHoles + + +def plotField (points,shape,title = "BoreField",SaveOrShow = True,Pdf = None): + ''' + Function graphs given bore hole field + Parameters + ------------- + :param points: [[float]] + array containing series of x,y points + :param shape: object shape + shape object representing obstruction used to generate series of x,y points + :param title: string + tile of figure + :param SaveOrShow: boolean + True if diplay is desired, False if figure should be saved + :param Pdf: Pdf Object + Object where figure should be saved + :return: Nothing + + ''' + plt.clf() + fig = plt.figure(figsize=(10,7)) + plt.rcParams.update({"font.size": 18}) + plt.title(title) + for i in range(len(shape.c)): + if i == len(shape.c)-1: + plt.plot([shape.c[i][0],shape.c[0][0]],[shape.c[i][1],shape.c[0][1]], "r") + else: + plt.plot([shape.c[i][0], shape.c[i+1][0]], [shape.c[i][1],shape.c[i+1][1]], "r") + X = points[:,0] + Y = points[:,1] + plt.plot(X,Y,"bo") + plt.axis("equal") + if SaveOrShow: + plt.show() + else: + Pdf.savefig(fig) + plt.close() + return + +def main(): + return + +def ProcessRows (BaseSLoc,rowspace,shape,XSpac,s1,BaseWidth,s2,i): + ''' + Function generates a row of the borefield + *Note: the formatting from the rows can be a little unexpected. Some adjustment + may be required to correct the formatting. The genBoreHoleConfig function already accounts for this. + Parameters + ------------- + :param BaseSLoc: [float,float] + x,y location of the leftmost point of the trapezoid + :param rowspace: float + the spacing between rows + :param shape: shape object + object representing "no-go" zone + :param XSpac: float + minimum spacing between columns + :param s1: float + slope of the left side of the trapezoid + :param BaseWidth: float + with of the base of the trapezoid + :param s2: float + slope of the rights side of the trapezoid + :param i: int + the index of the currunt row + :return: [[float]] + two dimensional array containing the x,y values of the bore holes for this row + ''' + rA = [] + currentY = BaseSLoc[1] + rowspace * i + rowsx = BaseSLoc[0] + s1 * currentY + rowex = BaseSLoc[0] + BaseWidth + s2 * currentY + currentX = rowsx + ncol = int((rowex - rowsx) // XSpac) + inters = shape.lineintersect([currentX, currentY, rowex, currentY]) + noin = len(inters) + if noin > 0: + inters.sort(key=operator.itemgetter(0)) + if noin > 1: + if inters[0][0] < rowsx and inters[noin - 1][0] > rowex: + inside = False + for inter in inters: + if inter[0] >= rowsx and inter[0] <= rowex: + inside = True + if not inside: + if shape.pointintersect([(rowex+rowsx)/2,currentY]): + return + inters = np.array(inters) + indices = [j for j in range(noin) if not (inters[j][0] > rowex or inters[j][0] < rowsx)] + inters = inters[indices] + noin = len(inters) + for i in range(noin-1): + spac = float(inters[i+1][0]-inters[i][0]) + if spac < XSpac and shape.pointintersect([(inters[i+1][0]+inters[i][0])/2,currentY]): + inters[i+1][0] += (XSpac-spac)/2 + inters[i][0] -= (XSpac-spac)/2 + if ncol < 1: + if not shape.pointintersect([(rowsx + rowex) / 2, currentY]): + return[[(rowsx + rowex) / 2, currentY]] + else: + if noin == 0: + distribute(rowsx, rowex, currentY, XSpac,rA) + elif noin == 2: + if shape.pointintersect([(inters[0][0] + inters[1][0]) / 2, currentY]): + distribute(rowsx, inters[0][0] , currentY, XSpac,rA) + distribute(inters[1][0], rowex, currentY, XSpac,rA) + else: + distribute(inters[0][0], inters[1][0], currentY, XSpac,rA) + elif noin == 1: + if shape.pointintersect([rowex, currentY]): + distribute(rowsx, inters[0][0], currentY, XSpac,rA) + else: + distribute(inters[0][0], rowex, currentY, XSpac,rA) + elif noin % 2 == 0: + i = 0 + while i < noin: + if i == 0: + distribute(rowsx, inters[0][0], currentY, XSpac,rA) + i = 1 + continue + elif i == noin - 1: + distribute(inters[noin - 1][0], rowex, currentY, XSpac,rA) + else: + distribute(inters[i][0], inters[i + 1][0], currentY, XSpac,rA) + i += 2 + else: + if shape.pointintersect([rowex, currentY]): + i = 0 + while i < noin: + if i == 0: + distribute(rowsx, inters[0][0], currentY, XSpac,rA) + i = 1 + continue + elif i == noin - 1: + i += 2 + continue + else: + distribute(inters[i][0], inters[i + 1][0], currentY, XSpac,rA) + i += 2 + else: + i = 0 + while i < noin: + if i == noin - 1: + distribute(inters[noin - 1][0], rowex, currentY, XSpac,rA) + else: + distribute(inters[i][0], inters[i + 1][0], currentY, XSpac,rA) + i += 2 + return rA +def Between(x1,x2,x3): + ''' + Function determines whenther x1 lies between x2 and x3 + Parameters + ------------- + :param x1: [float,float] + point 1 + :param x2: [float,float] + point 2 + :param x3: [float,float] + point 3 + :return: + ''' + if x1 >= x2 and x1 <= x3: + return True + return False +def distribute(x1,x2,y,spacing,r): + ''' + Function generates a series of boreholes between x1 and x2 + Parameters + ------------- + :param x1: float + left x value + :param x2: float + right x value + :param y: float + y value of row + :param spacing: float + spacing between columns + :param r: [[float]] + existing array of points + :return: + ''' + if x2-x1 < spacing: + r.append([(x2+x1)/2, y]) + return + currentX = x1 + actncol = int((x2 - x1) // spacing) + actSpac = float(x2 - x1) / actncol + while (currentX - x2) <= (1e-12): + r.append([currentX, y]) + currentX += actSpac + return + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/nogo.py b/nogo.py new file mode 100644 index 0000000..4b8a547 --- /dev/null +++ b/nogo.py @@ -0,0 +1,271 @@ +''' +*nogo.py* +creates and handles rectangular nogo zones + +Classes: + rec +Functions: + + angleBetween(float,float,float) - > float + vectorintersect()->[[float]] + +''' +from math import cos +from math import sin +from math import acos +from math import sqrt +from math import pi +import numpy as np +from rectangular_nogo import rec +import matplotlib.pyplot as plt +class Shapes: + ''' + a class to represent nogo zones + + Attributes + ---------- + cx : float + the x value of the centroid + cy : float + the y value of the centroid + xw: float + the xwidth + yw : float + the ywidth + theta: float + the rotation of the shape + c00: [float,float] + x,y location of 1 vertex + c01: [float,float] + x,y location of 2nd vertex + c10: [float,float] + x,y location of 3rd vertex + c11: [float,float] + x,y location of 4th vertex + c: [[float]] + array containing x,y locations of vertices + maxy : float + maximum y value of shape + miny: float + minimum y value of shape + maxx: float + maximum x value of shape + minx: float + minimum x value of shape + Methods + ------- + lineintersect(xy) + determines the intersection of this shape and the given line segment + pointintersect(xy) + determines whether the given point is inside of the rectangle + ''' + def __init__(self, cx, cy, xw, yw, theta=0,sh = "B"): + ''' + contructs a shape object + + Parameters + ---------- + :param cx: float + the x location of the centroid + :param cy: float + the y location of the centroid + :param xw: float + the width in the x dir + :param yw: float + the width in the y dir + :param theta: float + the amount of rotation in radians + :param sh: string + string specifying the desired shape supports: + B,S,L,U,T,BL + ''' + self.xw = xw + self.yw = yw + self.c = [] + self.rects = None + dy1 = .5*(yw+xw)*sin(pi/2 - theta) + dy2 = .5*(yw-xw)*sin(theta) + dx1 = .5*(yw+xw)*cos(pi/2-theta) + dx2 = .5*(yw-xw)*cos(theta) + if sh == "B": + self.rects = [rec(cx,cy,xw,yw,theta)] + self.c=self.rects[0].c + elif sh == "S": + self.rects = [rec(cx,cy,xw,yw,theta), + rec(cx-dx1+dx2,cy+dy1+dy2,yw,xw,theta), + rec(cx+dx1-dx2,cy-dy1-dy2,yw,xw,theta)] + self.c = [self.rects[2].c[0],self.rects[2].c[1],self.rects[0].c[2],self.rects[1].c[1],self.rects[1].c[2],self.rects[1].c[3],self.rects[0].c[0],self.rects[2].c[3]] + elif sh == "L": + self.rects = [rec(cx, cy, xw, yw, theta), + rec(cx + dx1+dx2, cy - dy1+dy2, yw, xw, + theta)] + self.c = [self.rects[1].c[0], self.rects[1].c[1], self.rects[1].c[2], self.rects[0].c[1], + self.rects[0].c[2], self.rects[0].c[3]] + elif sh == "U": + self.rects = [rec(cx, cy, xw, yw, theta), + rec(cx - dx1+dx2, + cy+dy1+dy2, yw, xw, theta), + rec(cx + dx1 + dx2, + cy - dy1 + dy2, yw, xw, theta)] + self.c = [self.rects[2].c[0], self.rects[2].c[1], self.rects[2].c[2], self.rects[0].c[1], self.rects[0].c[2], + self.rects[1].c[1], self.rects[1].c[2], self.rects[1].c[3]] + elif sh == "T": + self.rects = [rec(cx,cy,xw,yw,theta),rec(cx-dx1, cy+dy1, yw, xw, theta)] + self.c = [self.rects[0].c[0], self.rects[0].c[1], self.rects[0].c[2], self.rects[1].c[1], + self.rects[1].c[2], self.rects[1].c[3],self.rects[1].c[0],self.rects[0].c[3]] + + else: + self.rects = [rec(cx, cy, xw, yw, theta), + rec(cx + dx1 - dx2, + cy - dy1 - dy2, yw, xw, + theta)] + self.c = [self.rects[1].c[0], self.rects[1].c[1], self.rects[0].c[2], self.rects[0].c[3], + self.rects[0].c[0], self.rects[1].c[3]] + self.maxx = max(self.c[:][0]) + self.minx = min(self.c[:][0]) + self.maxy = max(self.c[:][1]) + self.miny = min(self.c[:][1]) + + def lineintersect(self, xy): + ''' + returns the intersections between a line segment and the rectanlge + + Parameters + ----------- + :param xy: [float,float,float,float] + the x,y values of both endpoints of the line segment + :return: [[float]] + the x,y values of the intersections + ''' + x1, y1, x2, y2 = xy + rA = [] + for i in range(len(self.c)): + if i == len(self.c)-1: + c1 = self.c[len(self.c)-1] + c2 = self.c[0] + r = vectorintersect([c1[0], c1[1], c2[0], c2[1]], [x1-100, y1, x2+100, y2]) + if r != None: + rA.append(r) + else: + c1 = self.c[i] + c2 = self.c[i + 1] + r = vectorintersect([c1[0], c1[1], c2[0], c2[1]], [x1-100, y1, x2+100, y2]) + if r != None: + rA.append(r) + return rA + + def pointintersect(self, xy): + ''' + returns whether the given point is inside of the rectangle + + Parameters + ----------- + :param xy: [float,float] + x,y value of point + :return: boolean + true if inside, false if not + ''' + for rect in self.rects: + if rect.pointintersect(xy): + return True + return False + def getArea (self): + ''' + returns area of shape + :return: float + area of shape + ''' + return len(self.rects)*self.xw*self.yw + + + +def angleBetween(a, b, c): + ''' + gives the angle opposite of side a + + Parameters + ----------- + :param a: float + length of side a + :param b: float + length of side b + :param c: float + length of side c + :return: double + angle across from a in triangle + ''' + if a == 0 or b == 0 or c == 0: + return 0 + r = (a * a + b * b - c * c) / (2 * a * b) + if r > 1 or r < -1: + return 0 + else: + return acos(r) + + +def vectorintersect(l1, l2): + ''' + gives the intersection between two line segments + + Parameters + ----------- + :param l1: [[float]] + endpoints of first line segment + :param l2: [[float]] + endpoints of the second line segment + :return: [float,float] + x,y values of intersection (returns None if there is none) + ''' + x11, y11, x12, y12 = l1 + x21, y21, x22, y22 = l2 + + p = np.array([x11, y11]) + q = np.array([x21, y21]) + r = np.array([x12 - x11, y12 - y11]) + s = np.array([x22 - x21, y22 - y21]) + rcs = np.cross(r, s) + if rcs != 0: + qpr = np.cross(np.subtract(q, p), r) + qps = np.cross(np.subtract(q, p), s) + u = qpr / rcs + t = qps / rcs + if 0 <= t <= 1 and 0 <= u <= 1: + return [p[0] + t * r[0], p[1] + t * r[1]] + return None +#Testing method +def main(): + ''' + tests shape class + + Parameters + ----------- + :return: none + ''' + rect = Shapes(3.0,4.0,4.0,9.0,pi/5.0,"S") + xy = [-3,11,10,11] + points = rect.lineintersect(xy) + for i in range(len(rect.c)): + if i == len(rect.c)-1: + plt.plot([rect.c[i][0],rect.c[0][0]],[rect.c[i][1],rect.c[0][1]], "r") + else: + plt.plot([rect.c[i][0], rect.c[i+1][0]], [rect.c[i][1], rect.c[i+1][1]], "r") + plt.plot([xy[0],xy[2]],[xy[1],xy[3]],"b") + plt.plot([-10,10],[rect.maxy,rect.maxy],"y") + X = np.zeros(len(points)) + Y = np.zeros(len(points)) + for i in range(len(points)): + X[i] = points[i][0] + Y[i] = points[i][1] + plt.plot(X,Y,"go") + plt.axis("equal") + plt.show() + for i in range(len(rect.rects)): + plt.plot([rect.rects[i].c00[0], rect.rects[i].c01[0]], [rect.rects[i].c00[1], rect.rects[i].c01[1]], "r") + plt.plot([rect.rects[i].c01[0], rect.rects[i].c11[0]], [rect.rects[i].c01[1], rect.rects[i].c11[1]], "r") + plt.plot([rect.rects[i].c11[0], rect.rects[i].c10[0]], [rect.rects[i].c11[1], rect.rects[i].c10[1]], "r") + plt.plot([rect.rects[i].c10[0], rect.rects[i].c00[0]], [rect.rects[i].c10[1], rect.rects[i].c00[1]], "r") + + plt.axis("equal") + plt.show() +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/rectangular_nogo.py b/rectangular_nogo.py new file mode 100644 index 0000000..634da60 --- /dev/null +++ b/rectangular_nogo.py @@ -0,0 +1,247 @@ +''' +*rectangular_nogo.py* +creates and handles rectangular nogo zones + +Classes: + rec +Functions: + + angleBetween(float,float,float) - > float + vectorintersect()->[[float]] + +''' +from math import cos +from math import sin +from math import acos +from math import sqrt +from math import pi +import numpy as np +import matplotlib.pyplot as plt + + +class rec: + ''' + a class to represent rectangular nogo zones + + Attributes + ---------- + cx : float + the x value of the centroid + cy : float + the y value of the centroid + xw: float + the xwidth + yw : float + the ywidth + theta: float + the rotation of the rectangle + c00: [float,float] + x,y location of 1 vertex + c01: [float,float] + x,y location of 2nd vertex + c10: [float,float] + x,y location of 3rd vertex + c11: [float,float] + x,y location of 4th vertex + c: [[float]] + array containing x,y locations of vertices + maxy : float + maximum y value of shape + miny: float + minimum y value of shape + maxx: float + maximum x value of shape + minx: float + minimum x value of shape + Methods + ------- + lineintersect(xy) + determines the intersection of this shape and the given line segment + pointintersect(xy) + determines whether the given point is inside of the rectangle + ''' + def __init__(self,cx,cy,xw,yw,theta=0): + ''' + contructs a rectangle object + + Parameters + ---------- + :param cx: float + the x location of the centroid + :param cy: float + the y location of the centroid + :param xw: float + the width in the x dir + :param yw: float + the width in the y dir + :param theta: float + the amount of rotation in radians + ''' + self.cx = cx + self.cy = cy + self.xw = xw + self.yw = yw + self.theta = theta + hxwc = .5*xw*cos(theta) + hywc = .5*yw*cos(theta) + hxws = .5*xw*sin(theta) + hyws = .5*yw*sin(theta) + self.c00 =[ cx - hxwc + hyws,cy - hxws-hywc] + self.c10 =[ cx +hxwc+hyws,cy + hxws-hywc] + self.c01 =[ cx -hxwc-hyws,cy - hxws+hywc] + self.c11 =[ cx + hxwc - hyws,cy + hxws + hywc] + self.c = [self.c00,self.c10,self.c11,self.c01] + self.maxy = max(self.c00[1],self.c10[1],self.c01[1],self.c11[1]) + self.miny = min(self.c00[1],self.c10[1],self.c01[1],self.c11[1]) + self.maxx = max(self.c00[0],self.c10[0],self.c01[0],self.c11[0]) + self.minx = min(self.c00[0],self.c10[0],self.c01[0],self.c11[0]) + def lineintersect(self,xy): + ''' + returns the intersections between a line segment and the rectanlge + + Parameters + ----------- + :param xy: [float,float,float,float] + the x,y values of both endpoints of the line segment + :return: [[float]] + the x,y values of the intersections + ''' + x1,y1,x2,y2 = xy + if (y1 > self.maxy and y2 > self.maxy) or (y1 < self.miny and y2 < self.miny): + return [] + rA = [] + for i in range(4): + if i == 3: + c1 = self.c[3] + c2 = self.c[0] + r = vectorintersect([c1[0],c1[1],c2[0],c2[1]],[x1,y1,x2,y2]) + if r != None: + rA.append(rA, r,) + else: + c1 = self.c[i] + c2 = self.c[i+1] + r = vectorintersect([c1[0], c1[1], c2[0], c2[1]], [x1, y1, x2, y2]) + if r != None: + rA.append(rA, r,) + return rA + + + def pointintersect(self,xy): + ''' + returns whether the given point is inside of the rectangle + + Parameters + ----------- + :param xy: [float,float] + x,y value of point + :return: boolean + true if inside, false if not + ''' + x,y = xy + if y > self.maxy or y 1 or r < -1: + return 0 + else: + return acos(r) + + +def vectorintersect(l1,l2): + ''' + gives the intersection between two line segments + + Parameters + ----------- + :param l1: [[float]] + endpoints of first line segment + :param l2: [[float]] + endpoints of the second line segment + :return: [float,float] + x,y values of intersection (returns None if there is none) + ''' + x11,y11,x12,y12 = l1 + x21,y21,x22,y22 = l2 + + p = np.array([x11,y11]) + q = np.array([x21,y21]) + r = np.array([x12-x11,y12-y11]) + s = np.array([x22 - x21, y22 - y21]) + rcs = np.cross(r,s) + if rcs != 0: + qpr = np.cross(np.subtract(q, p), r) + qps = np.cross(np.subtract(q, p), s) + u = qpr/rcs + t=qps/rcs + if 0 <= t <= 1 and 0 <= u <= 1: + return [p[0]+t*r[0], p[1] + t*r[1]] + return None + +def main(): + ''' + tests rectangle class + + Parameters + ----------- + :return: none + ''' + rect = rec(3.0,4.0,7.0,9.0,pi/5.0) + print("Should be True: ",rect.pointintersect([4.0,5.0])) + print("Should be False: ", rect.pointintersect([-3.0, 9.0])) + print("Should be False: ", rect.pointintersect([4.0, -1.0])) + print("Should be True: ", rect.pointintersect([-2.0, 5.5])) + xy = [-3,-2,9,7] + points = rect.lineintersect(xy) + plt.plot([rect.c00[0], rect.c01[0]], [rect.c00[1], rect.c01[1]], "r") + plt.plot([rect.c01[0], rect.c11[0]], [rect.c01[1], rect.c11[1]], "r") + plt.plot([rect.c11[0], rect.c10[0]], [rect.c11[1], rect.c10[1]], "r") + plt.plot([rect.c10[0], rect.c00[0]], [rect.c10[1], rect.c00[1]], "r") + plt.plot([xy[0],xy[2]],[xy[1],xy[3]],"b") + X = np.zeros(len(points)) + Y = np.zeros(len(points)) + for i in range(len(points)): + X[i] = points[i][0] + Y[i] = points[i][1] + plt.plot(X,Y,"go") + plt.show() +if __name__ == "__main__": + main() \ No newline at end of file