I am using pythom 3.6 and I have came across a problem. There are given n points by x and y coordinates (n is from 4 to 200) and I need to find the 4 points from those n which form the biggest general quadrangle (any convex shape formed by 4 points).
I can think of a solution including 4 for cycles with calculating the area of the quadrangle given by point in the for cycles, but it is extremaly slow. Do you know about anything faster?
The point are given like this:
B = np.array([[ 1., 2.], [ 0., 0.], [ 1., 3.], [ -5., 6.], [ -6., 3.], [ 1., 5.], [ -1., 2.], [ 1., -3.], [ 4., 2.]])
The next level is when I get N points given by x, y and z coordinates (N is between 8 and 500) and I should find the biggest (in volume) hexahedron (the shape defined by 8 points) - I have no idea of the solution.
There is no need for right angles, just shapes defined by 4 (8) points. Any suggestions?
Background: I have quite complex 3D models of building which I need to simplify to one specific program for computations. The details on the building are not needed. All information about the buildings is in file.obj exported from Blender.
So, the answer for the 2D and 3D problem is similar. Because it are buildings, we can split the 3D model to two 2D models with the base of the building and the roof (and anything between is considered as roof too).
Then we are searching for the quadrangle to approximate the points in a (approximately) plane. We need to find the centre (not by mean from all points, but (max+min)/2 in both directons). We move the origin to the centre by calculation of the points - centre. The points then should be devided by quadrants (x>0 & y>0, x<0, & y>0, x<0, & y<0, x>0, & y<0) and for each quadrant we calculate the most distant point (if there is nan, we take the origin [0,0]).
Using the Shoelace formula we calculate the area taken by those 4 points from each quadrants keeping that value. After that we rotate the points around the origin by 1 degree up to 90 degrees. Each time we calculate the area and we are searching for the maximum area. The points shoving the maximum area are the desired points. The code (I know it is not smooth code, can be optimized. But it works!!):
def getCorners(points):
maxPoint = np.max(points[:,0])
mayPoint = np.max(points[:,1])
minPoint = np.min(points[:,0])
miyPoint = np.min(points[:,1])
meanPoint = np.array([(maxPoint + minPoint)/2, (mayPoint + miyPoint)/2])
normPoints = points[:,0:2] - meanPoint[0:2]
areaMaximum = -1
finID1 = 0
finID2 = 1
finID3 = 2
finID4 = 3
numrot = 360
for alpha in range(0,numrot):
topright = np.where((normPoints[:,0] > 0) & (normPoints[:,1] > 0))
topleft = np.where((normPoints[:,0] < 0) & (normPoints[:,1] > 0))
bottomleft = np.where((normPoints[:,0] < 0) & (normPoints[:,1] < 0))
bottomright = np.where((normPoints[:,0] > 0) & (normPoints[:,1] < 0))
q1 = normPoints[topright]
q2 = normPoints[topleft]
q3 = normPoints[bottomleft]
q4 = normPoints[bottomright]
if len(q1) == 0:
q1 = np.array([[0,0],[0,0]])
if len(q2) == 0:
q2 = np.array([[0,0],[0,0]])
if len(q3) == 0:
q3 = np.array([[0,0],[0,0]])
if len(q4) == 0:
q4 = np.array([[0,0],[0,0]])
D1 = q1[:,0]*q1[:,0] + q1[:,1]*q1[:,1]
D2 = q2[:,0]*q2[:,0] + q2[:,1]*q2[:,1]
D3 = q3[:,0]*q3[:,0] + q3[:,1]*q3[:,1]
D4 = q4[:,0]*q4[:,0] + q4[:,1]*q4[:,1]
ID1 = np.argmax(D1)
ID2 = np.argmax(D2)
ID3 = np.argmax(D3)
ID4 = np.argmax(D4)
vertices = [[q1[ID1,0],q1[ID1,1]],[q2[ID2,0],q2[ID2,1]],[q3[ID3,0],q3[ID3,1]],[q4[ID4,0],q4[ID4,1]]]
area = polygonArea(vertices)
if area > areaMaximum:
areaMaximum = area
if len(topright[0]) == 0:
finID1 = 0
else:
finID1 = topright[0][ID1]
if len(topleft[0]) == 0:
finID2 = 0
else:
finID2 = topleft[0][ID2]
if len(bottomleft[0]) == 0:
finID3 = 0
else:
finID3 = bottomleft[0][ID3]
if len(bottomright[0]) == 0:
finID4 = 0
else:
finID4 = bottomright[0][ID4]
# rotate
for opi in range(0,len(normPoints)):
normPoints[opi] = rotate_origin_only(normPoints[opi], 90/numrot/180*np.pi)
return [finID1,finID2,finID3,finID4]
def rotate_origin_only(xy, radians):
"""Only rotate a point around the origin (0, 0)."""
x, y = xy
xx = x * math.cos(radians) + y * math.sin(radians)
yy = -x * math.sin(radians) + y * math.cos(radians)
return xx, yy
def polygonArea(vertices):
#A function to apply the Shoelace algorithm
numberOfVertices = len(vertices)
sum1 = 0
sum2 = 0
for i in range(0,numberOfVertices-1):
sum1 = sum1 + vertices[i][0] * vertices[i+1][1]
sum2 = sum2 + vertices[i][1] * vertices[i+1][0]
#Add xn.y1
sum1 = sum1 + vertices[numberOfVertices-1][0]*vertices[0][1]
#Add x1.yn
sum2 = sum2 + vertices[0][0]*vertices[numberOfVertices-1][1]
area = abs(sum1 - sum2) / 2
return area