【发布时间】：2022-01-21 15:13:18
【问题描述】：
我正在尝试修改 maxrect 库以解决不受方向限制的最大矩形。
查看代码我看到的约束是：
""" :param coordinates:
A list of of [x, y] pairs describing a closed, convex polygon.
"""
coordinates = np.array(coordinates)
x_range = np.max(coordinates, axis=0)[0]-np.min(coordinates, axis=0)[0]
y_range = np.max(coordinates, axis=0)[1]-np.min(coordinates, axis=0)[1]
scale = np.array([x_range, y_range])
sc_coordinates = coordinates/scale
poly = Polygon(sc_coordinates)
inside_pt = (poly.representative_point().x,
poly.representative_point().y)
A1, A2, B = pts_to_leq(sc_coordinates)
bl = cvxpy.Variable(2)
tr = cvxpy.Variable(2)
br = cvxpy.Variable(2)
tl = cvxpy.Variable(2)
obj = cvxpy.Maximize(cvxpy.log(tr[0] - bl[0]) + cvxpy.log(tr[1] - bl[1]))
constraints = [bl[0] == tl[0],
br[0] == tr[0],
tl[1] == tr[1],
bl[1] == br[1],
]
for i in range(len(B)):
if inside_pt[0] * A1[i] + inside_pt[1] * A2[i] <= B[i]:
constraints.append(bl[0] * A1[i] + bl[1] * A2[i] <= B[i])
constraints.append(tr[0] * A1[i] + tr[1] * A2[i] <= B[i])
constraints.append(br[0] * A1[i] + br[1] * A2[i] <= B[i])
constraints.append(tl[0] * A1[i] + tl[1] * A2[i] <= B[i])
else:
constraints.append(bl[0] * A1[i] + bl[1] * A2[i] >= B[i])
constraints.append(tr[0] * A1[i] + tr[1] * A2[i] >= B[i])
constraints.append(br[0] * A1[i] + br[1] * A2[i] >= B[i])
constraints.append(tl[0] * A1[i] + tl[1] * A2[i] >= B[i])
也就是说，他们将外接多边形中的每个点转换为 Ax + Ay = B 并检查矩形的角是否在其中，并最大化对角线。此外，还有 4 个约束可确保角的角度与参考框架对齐。
我在想我可以删除这 4 个约束。
但是，这允许矩形超出外接多边形的边界。这可能意味着我对上述约束的目的并不完全正确。也可能是现在允许这些点在参考框架内偏离它们的基数方向，所以我尝试添加不同的约束来保持矩形点的基数：
if aligned:
constraints.append(bl[0] == tl[0])
constraints.append(br[0] == tr[0])
constraints.append(tl[1] == tr[1])
constraints.append(bl[1] == br[1])
else:
constraints.append(bl[0] < br[0])
constraints.append(bl[0] < tr[0])
constraints.append(tl[0] < br[0])
constraints.append(tl[0] < tr[0])
constraints.append(bl[1] < tl[1])
constraints.append(bl[1] < tr[1])
constraints.append(br[1] < tl[1])
constraints.append(br[1] < tr[1])
有人可以帮我看看我缺少什么吗？
示例问题：
square = Polygon([(0,0), (0,2), (2,2), (2,0)], [
[(0.5,0.5), (0.5,1.2), (1,1.2), (1,0.5)],
[(1.1,1.1), (1.1,1.5), (1.5,1.5), (1.5,1.1)]
])
line = LineString((square.exterior.coords[0], square.exterior.coords[1]))
max_hull = find_maximal_convex_hull(line, square)
print('max hull', max_hull.wkt, max_hull.area)
# this line calls my cvypy script
pa = get_maximal_rectangle(max_hull.exterior.coords)
max_rect = Polygon([(pa[0][0],pa[0][1]), (pa[0][0],pa[1][1]), (pa[1][0],pa[1][1]), (pa[1][0],pa[0][1])])
print('max rectangle', pa)
plt.axis(xmin=-0.125,xmax=2.125,ymin=-0.125,ymax=2.125)
plt.plot(*square.exterior.xy, color='g')
[plt.plot(*i.xy, color='y') for i in square.interiors]
plt.plot(*max_hull.exterior.xy, color='b')
[plt.plot(*i.xy, color='r') for i in max_hull.interiors]
plt.plot(*max_rect.exterior.xy, color='r')
plt.show()
在上面的代码中，我得到了一个复杂的形状和一条边。我必须将其切割成（大致/几乎）最大的凸形，并且边缘必须与该形状相交。这项工作我已经完成了，它被分配给变量max_hull。
现在，无论方向如何，我都想要形状中最大的矩形。在下图中，根据上面的代码绘制，我用绿色显示了外部正方形（左边缘是所需的边缘），它的孔用黄色显示，剩余的凸形，即给我的 cvypy 的外部多边形用蓝色显示， -- 它从可见的底线一直延伸到顶部 -- 并且从 cvxpy 返回的候选矩形为红色。
日志输出为：
<ipython-input-2-efd2419ff139>:272: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
for base in split(S, ll):
455 unique polygons considered
max hull POLYGON ((0 0.95, 0 2, 2 2, 2 1.95, 0 0.95)) 1.0999999999999999
max rectangle (array([7.13062842e-09, 9.50000011e-01]), array([1.99999998, 1.99999999]))
为了在现代 cvypy 和 python 中使用 maxrect 库，请将 get_max_rectangle（在 init.py 中）solve 语句更改为 prob.solve() 和 return 语句。应用我的代码后，这是我正在使用的自定义函数：
import numpy as np
import cvxpy
from shapely.geometry import Polygon
def rect2poly(ll, ur):
"""
Convert rectangle defined by lower left/upper right
to a closed polygon representation.
"""
x0, y0 = ll
x1, y1 = ur
return [
[x0, y0],
[x0, y1],
[x1, y1],
[x1, y0],
[x0, y0]
]
def get_intersection(coords):
"""Given an input list of coordinates, find the intersection
section of corner coordinates. Returns geojson of the
interesection polygon.
"""
ipoly = None
for coord in coords:
if ipoly is None:
ipoly = Polygon(coord)
else:
tmp = Polygon(coord)
ipoly = ipoly.intersection(tmp)
# close the polygon loop by adding the first coordinate again
first_x = ipoly.exterior.coords.xy[0][0]
first_y = ipoly.exterior.coords.xy[1][0]
ipoly.exterior.coords.xy[0].append(first_x)
ipoly.exterior.coords.xy[1].append(first_y)
inter_coords = zip(
ipoly.exterior.coords.xy[0], ipoly.exterior.coords.xy[1])
inter_gj = {"geometry":
{"coordinates": [inter_coords],
"type": "Polygon"},
"properties": {}, "type": "Feature"}
return inter_gj, inter_coords
def two_pts_to_line(pt1, pt2):
"""
Create a line from two points in form of
a1(x) + a2(y) = b
"""
pt1 = [float(p) for p in pt1]
pt2 = [float(p) for p in pt2]
try:
slp = (pt2[1] - pt1[1]) / (pt2[0] - pt1[0])
except ZeroDivisionError:
slp = 1e5 * (pt2[1] - pt1[1])
a1 = -slp
a2 = 1.
b = -slp * pt1[0] + pt1[1]
return a1, a2, b
def pts_to_leq(coords):
"""
Converts a set of points to form Ax = b, but since
x is of length 2 this is like A1(x1) + A2(x2) = B.
returns A1, A2, B
"""
A1 = []
A2 = []
B = []
for i in range(len(coords) - 1):
pt1 = coords[i]
pt2 = coords[i + 1]
a1, a2, b = two_pts_to_line(pt1, pt2)
A1.append(a1)
A2.append(a2)
B.append(b)
return A1, A2, B
def get_maximal_rectangle(coordinates, aligned=False):
"""
Find the largest, inscribed, axis-aligned rectangle.
:param coordinates:
A list of of [x, y] pairs describing a closed, convex polygon.
"""
coordinates = np.array(coordinates)
x_range = np.max(coordinates, axis=0)[0]-np.min(coordinates, axis=0)[0]
y_range = np.max(coordinates, axis=0)[1]-np.min(coordinates, axis=0)[1]
scale = np.array([x_range, y_range])
sc_coordinates = coordinates/scale
rep_point = Polygon(sc_coordinates).representative_point()
inside_pt = (rep_point.x, rep_point.y)
A1, A2, B = pts_to_leq(sc_coordinates)
bl = cvxpy.Variable(2)
tr = cvxpy.Variable(2)
br = cvxpy.Variable(2)
tl = cvxpy.Variable(2)
obj = cvxpy.Maximize(cvxpy.log(tr[0] - bl[0]) + cvxpy.log(tr[1] - bl[1]))
constraints = []
if aligned:
constraints.append(bl[0] == tl[0])
constraints.append(br[0] == tr[0])
constraints.append(tl[1] == tr[1])
constraints.append(bl[1] == br[1])
else:
constraints.append(bl[0] < br[0])
constraints.append(bl[0] < tr[0])
constraints.append(tl[0] < br[0])
constraints.append(tl[0] < tr[0])
constraints.append(bl[1] < tl[1])
constraints.append(bl[1] < tr[1])
constraints.append(br[1] < tl[1])
constraints.append(br[1] < tr[1])
for i in range(len(B)):
if inside_pt[0] * A1[i] + inside_pt[1] * A2[i] <= B[i]:
constraints.append(bl[0] * A1[i] + bl[1] * A2[i] <= B[i])
constraints.append(tr[0] * A1[i] + tr[1] * A2[i] <= B[i])
constraints.append(br[0] * A1[i] + br[1] * A2[i] <= B[i])
constraints.append(tl[0] * A1[i] + tl[1] * A2[i] <= B[i])
else:
constraints.append(bl[0] * A1[i] + bl[1] * A2[i] >= B[i])
constraints.append(tr[0] * A1[i] + tr[1] * A2[i] >= B[i])
constraints.append(br[0] * A1[i] + br[1] * A2[i] >= B[i])
constraints.append(tl[0] * A1[i] + tl[1] * A2[i] >= B[i])
prob = cvxpy.Problem(obj, constraints)
#prob.solve(solver=cvxpy.CVXOPT, verbose=False, max_iters=1000, reltol=1e-9)
#prob.solve(solver=cvxpy.SCS, verbose=True, use_indirect=False, max_iters=int(1e5))
prob.solve()
bottom_left = np.array(bl.value).T * scale
top_right = np.array(tr.value).T * scale
#return list(bottom_left[0]), list(top_right[0])
return bottom_left, top_right
上面没有修改其他功能。
【问题讨论】：