有没有快速的公式可以找到由N个点给出的一般四边形(六面体)的最大面积?

Is there any fast formula to find the biggest area of a general quadrangle (hexahedron) given by N points?

我正在使用 pythom 3.6,但遇到了一个问题。给定了 x 和 y 坐标的 n 个点(n 是从 4 到 200),我需要从这 n 个点中找到形成最大一般四边形(由 4 个点形成的任何凸形)的 4 个点。

我可以想到一个解决方案,包括 4 个 for 循环,计算 for 循环中点给定的四边形的面积,但它非常慢。你知道什么更快吗?

积分是这样的:

B = np.array([[ 1., 2.], [ 0., 0.], [ 1., 3.], [ -5., 6.], [ -6., 3.], [ 1., 5.], [ -1., 2.], [ 1., -3.], [ 4., 2.]])

下一级是当我得到由 x、y 和 z 坐标给出的 N 个点(N 在 8 到 500 之间)时,我应该找到最大的(体积)六面体(由 8 个点定义的形状)-我不知道解决方案。

不需要直角,只需由 4 (8) 个点定义的形状。有什么建议吗?


背景: 我有相当复杂的 3D 建筑模型,我需要将其简化为一个特定的计算程序。不需要有关建筑物的详细信息。关于建筑物的所有信息都在从 Blender 导出的 file.obj 中。

构建 convex hull 个所有点。

然后求顶点属于壳的最大面积四边形。 如果船体数 N 很小,您可以只检查所有对角线。

否则考虑使用更高级的算法,例如:Maximum-Area Quadrilateral in a Convex Polygon, Revisited

所以,2D 和 3D 问题的答案是相似的。因为它是建筑物,所以我们可以将 3D 模型拆分为两个 2D 模型,包括建筑物的底部和屋顶(任何介于两者之间的都被视为屋顶)。

然后我们正在搜索四边形来近似(近似)平面中的点。我们需要找到中心(不是通过所有点的平均值,而是在两个方向上的 (max+min)/2)。我们通过计算点 - 中心将原点移动到中心。然后,这些点应按象限 (x>0 & y>0, x<0, & y>0, x<0, & y<0, x>0, & y<0) 划分,对于每个象限,我们计算最远的点(如果有 nan,我们取原点 [0,0])。

我们使用 Shoelace 公式计算每个象限中这 4 个点所占的面积,并保持该值。 之后,我们将点围绕原点旋转 1 度,最多旋转 90 度。每次我们计算面积,我们都在寻找最大面积。推动最大面积的点是所需的点。 代码(我知道这不是流畅的代码,可以优化。但它有效!!):

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