python - cv.MinAreaRect2 和 ArcGIS(GIS 软件)之间的结果差异。可能的错误？

Question

我有一组从多边形派生的点。我正在测试几种解决方案以获得最小面积或矩形。作为基准，我使用 ArcGIS (10.1)。

points = [(560036.4495758876, 6362071.890493258),
          (560036.4495758876, 6362070.890493258),
          (560036.9495758876, 6362070.890493258),
          (560036.9495758876, 6362070.390493258),
          (560037.4495758876, 6362070.390493258),
          (560037.4495758876, 6362064.890493258),
          (560036.4495758876, 6362064.890493258),
          (560036.4495758876, 6362063.390493258),
          (560035.4495758876, 6362063.390493258),
          (560035.4495758876, 6362062.390493258),
          (560034.9495758876, 6362062.390493258),
          (560034.9495758876, 6362061.390493258),
          (560032.9495758876, 6362061.390493258),
          (560032.9495758876, 6362061.890493258),
          (560030.4495758876, 6362061.890493258),
          (560030.4495758876, 6362061.390493258),
          (560029.9495758876, 6362061.390493258),
          (560029.9495758876, 6362060.390493258),
          (560029.4495758876, 6362060.390493258),
          (560029.4495758876, 6362059.890493258),
          (560028.9495758876, 6362059.890493258),
          (560028.9495758876, 6362059.390493258),
          (560028.4495758876, 6362059.390493258),
          (560028.4495758876, 6362058.890493258),
          (560027.4495758876, 6362058.890493258),
          (560027.4495758876, 6362058.390493258),
          (560026.9495758876, 6362058.390493258),
          (560026.9495758876, 6362057.890493258),
          (560025.4495758876, 6362057.890493258),
          (560025.4495758876, 6362057.390493258),
          (560023.4495758876, 6362057.390493258),
          (560023.4495758876, 6362060.390493258),
          (560023.9495758876, 6362060.390493258),
          (560023.9495758876, 6362061.890493258),
          (560024.4495758876, 6362061.890493258),
          (560024.4495758876, 6362063.390493258),
          (560024.9495758876, 6362063.390493258),
          (560024.9495758876, 6362064.390493258),
          (560025.4495758876, 6362064.390493258),
          (560025.4495758876, 6362065.390493258),
          (560025.9495758876, 6362065.390493258),
          (560025.9495758876, 6362065.890493258),
          (560026.4495758876, 6362065.890493258),
          (560026.4495758876, 6362066.890493258),
          (560026.9495758876, 6362066.890493258),
          (560026.9495758876, 6362068.390493258),
          (560027.4495758876, 6362068.390493258),
          (560027.4495758876, 6362068.890493258),
          (560027.9495758876, 6362068.890493258),
          (560027.9495758876, 6362069.390493258),
          (560028.4495758876, 6362069.390493258),
          (560028.4495758876, 6362069.890493258),
          (560033.4495758876, 6362069.890493258),
          (560033.4495758876, 6362070.390493258),
          (560033.9495758876, 6362070.390493258),
          (560033.9495758876, 6362070.890493258),
          (560034.4495758876, 6362070.890493258),
          (560034.4495758876, 6362071.390493258),
          (560034.9495758876, 6362071.390493258),
          (560034.9495758876, 6362071.890493258),
          (560036.4495758876, 6362071.890493258)]

一种解决方案是使用 OpenCV 的 cv.MinAreaRect2()。

The function cv.MinAreaRect2 finds a circumscribed rectangle of the minimal area for 2D point set by building convex hull for the set and applying rotating calipers technique to the hull.

import cv
# (x,y) - center point of the box
# (w,h) - width and height of the box
# theta - angle of rotation
((x,y),(w,h),th) = cv.MinAreaRect2(points)
print ((x,y),(w,h),th)
((560029.3125, 6362065.5), (10.28591251373291, 18.335756301879883), -63.43495178222656)
# get vertex
box_vtx = cv.BoxPoints(((x,y),(w,h),th))
print box_vtx 
((560035.1875, 6362074.0), (560018.8125, 6362066.0), (560023.4375, 6362057.0), (560039.8125, 6362065.0)

当我在 shapefile 中转换 box_vtx 以便在 ArcGIS 中查看并与 Minimum Bounding Geometry (Data Management) 进行比较时我可以看到这种差异，如下图所示， 其中:

• 红色 = 多边形边框
• 蓝色 = ArGIS 中面积最小的矩形 (10.1)
• 黄色和黑色 = OpenCV 中面积最小的矩形

使用 OpenCV 与此 post 中提出的解决方案进行比较:

import osgeo.gdal, ogr
import cv

poly = "...\\polygon.shp"
shp = osgeo.ogr.Open(poly)
layer = shp.GetLayer()
feature = layer.GetFeature(0)
geometry = feature.GetGeometryRef()
pts = geometry.GetGeometryRef(0)
# get point of the polygon border (the points above)
points = []
for p in xrange(pts.GetPointCount()):
    points.append((pts.GetX(p), pts.GetY(p)))
# Convex Hull
CH1 = geometry.ConvexHull
# i didn't find a method to extarct  the points
print CH1()
# works with openCV
cvxHull = cv.ConvexHull2(points, cv.CreateMemStorage(), return_points=True) 
print cvxHull
<cv2.cv.cvseq object at 0x00000000067CCF90>

Answer 1

我没有看过 OpenCV 代码，但在计算中似乎可能会从其他大 x,y 乘积中减去大 x,y 乘积。 x 值中的偏移量使用大约 19 位，y 值中的偏移量大约使用 23 位，因此这种减法可能会导致典型 double 携带的 53 位丢失大约 42 位。 (参见维基百科中的 floating point。)黄色和黑色矩形的大小和形状看起来合理，但显示的宽度和长度 (10.285..., 18.335...) 与 (10.393, 18.037) 存在 1% 的差异， related question 中出现的宽度和长度.

简而言之，OpenCV 在MinAreaRect2() 中可能存在舍入、下溢或溢出问题。或在 BoxPoints() .

需要检查或测试的内容包括:
• 显示并打印OpenCV凸包计算的点
• 在图表上，显示 OpenCV 和 ArcGIS 矩形的中心，并打印从这些中心到所显示矩形的角的距离
• 使用转换后的数据集重新运行计算和图表(见下文)，以减少彼此相减大数时的下溢效应

在 OpenCV 计算之前转换数据集可以将丢失的位数减少一半。通过如下代码生成一组新的点，并尝试使用新数据集进行计算:

s = points     # points = original data set
n = len(s)
cx = sum(zip(*s)[0])/n
cy = sum(zip(*s)[1])/n
points = map(lambda p: (p[0]-cx, p[1]-cy), s)
# Now points = translated data set

在上面的代码中，zip(*s)将 (x,y) 点列表解压缩为两个列表，其中 x 值在 zip(*s)[0] 中列出，y 值在 zip(*s)[1] 中列出。因此 (cx,cy) 表示 s 中列出的点的质心。映射表达式将函数应用于可迭代的元素。该函数返回由 (-cx,-cy) 转换的点。 map 表达式的值是平移点的列表。请注意，最好设置 cx = int(sum(zip(*s)[0])/n)和cy = int(sum(zip(*s)[1])/n)这样，如果您正在计算的话，格点仍然是格点。消除大偏移量的有益效果仍然存在，并且平移过程中会发生较少的舍入。

注意 - 我运行了从数据集和船体中减去偏移量的测试(假设船体如 question #13542855 或 question #13553884 所示)，并得到了不一致的结果与我在 #13542855 中给出的方法的结果不太一致。因此，很难根据您的数字确定问题发生的位置。下面显示了一个更简单的测试用例。该测试清楚地表明，较大的偏移会导致精度较差。也许您可以通过 ArcGIS 运行类似的测试数据。以下是测试结果的一半。 MinAreaRect2() 和 BoxPoints() 被赋予添加了偏移量的基数；对于显示的结果，计算后会减去偏移量，以便于比较结果。理想情况下，每列(除了 ox,oy 列)中的所有数字都相同；但随着 ox,oy 的增加，它们开始有所不同，并很快变得几乎是随机的。

      ox         oy      T cx      T cy    height     width     theta
       0          0    6.2500    7.7500   11.5709   13.7281  -78.6901
      64        125    6.2500    7.7500   11.5709   13.7281  -78.6901
     256        625    6.2500    7.7501   11.5709   13.7281  -78.6901
    1024       3125    6.2500    7.7500   11.5709   13.7281  -78.6901
    4096      15625    6.2500    7.7510   11.5709   13.7281  -78.6901
   16384      78125    6.2500    7.7578   11.5709   13.7281  -78.6901
   65536     390625    6.2500    7.7812   11.5709   13.7281  -78.6901
  262144    1953125    6.3125    7.7500   11.5709   13.7281  -78.6901
 1048576    9765625    6.2500    8.0000   11.5709   13.7281  -78.6901
 4194304   48828125    7.0000   11.0000   12.0000   14.0000  -90.0000
16777216  244140625    8.0000   15.0000   16.0000   14.0000  -90.0000
67108864 1220703125    8.0000  -21.0000   16.0000    0.0000   -0.0000

这是生成上面显示的数据的代码:

#!/usr/bin/python
import cv
base = [(1,1), (0,4), (2,9), (5,11), (8,14), (13,9), (14,4), (12,3), (2,1), (1,1)]
ox, oy, boxes = 0, 0, []
print '        ox         oy      T cx      T cy    height     width     theta'
for i in range(3,15):
    poly = map(lambda p: (p[0]+ox, p[1]+oy), base)
    (x,y), (w,h), th = cv.MinAreaRect2(poly)
    boxes.append((ox, oy, cv.BoxPoints(((x,y),(w,h),th))))
    print ('{:10d} {:10d} {:9.4f} {:9.4f} {:9.4f} {:9.4f} {:9.4f}'.format(
            ox, oy, x-ox, y-oy, w, h, th))
    ox, oy = 4**i, 5**i

print '\n        ox         oy       T x       T y       T x       T y       T x       T y       T x       T y'
for (ox, oy, box) in boxes:
    print '{:10d} {:10d}'.format(ox, oy),
    for p in box:
        print '{:9.4f} {:9.4f}'.format(p[0]-ox, p[1]-oy),
    print