cv.MinAreaRect2和ArcGIS（GIS软件）之间的差异的结果。 可能的错误？(Re

2019-08-03 05:30发布

站内文章 / 移动开发
154 0 0

我有一组从多边形派生点。 我为了获得最小的区域或矩形测试几个解决方案。 作为基准我使用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)]

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

功能cv.MinAreaRect2发现用于通过构建针对该组凸包并施加旋转卡钳技术在船体设置2D点的最小区域的外接矩形。 

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)

当我转换box_vtx在shape文件以查看在ArcGIS并用比较最小边界几何（数据管理）我可以看到如下面的图，其中示出这种差异：

  • 多边形的红色=边界
  • 蓝色=矩形与来自ArGIS最小区域（10.1）
  • 黄色和黑色=矩形与来自OpenCV的最小面积

与OpenCV的，在此提出的解决方案相比，工作职位 ： 

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，所以这样的减法可导致从53约42比特的损失典型的双精度数携带。 （参见浮点在维基百科。）的黄色和黑色矩形的大小和形状看起来合理的，但所显示的宽度和长度，（10.285 ...，18.335 ...）是从（10.393，18.037）1％不同，则宽度和长度中出现的一个相关的问题 。

总之，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值（* S）[0]和在拉链列出的y值（* s）实施1]。 因此（CX，CY）表示的上市点的质量中心。 地图表达应用一个函数可迭代的元素。 该函数返回由（-cx，-cy）翻译的点。 地图表达式的值是翻译的点的列表。 注意，这可能是最好设定cx = int(sum(zip(*s)[0])/n)和cy = int(sum(zip(*s)[1])/n)使得晶格点保持格点，如果那是你与计算。 去除大偏移的有益作用仍然存在，并且翻译过程中会发生更少的舍入。

注-我跑了从数据集中，然后从船体减去偏移测试（假设是如图所示的船体问题＃13542855 ，也许在问题＃135538​​84 ），并得到了不从我给方法的结果也一致不一致的结果在＃13542855。 因此，它是困难的问题出在哪里是基于你的号码出现的牵制。 一个更简单的测试案例如下所示。 这个试验清楚地表明，大偏移导致精度降低。 也许你可以通过ArcGIS的运行类似的测试数据。 以下是从测试结果的一半。 MinAreaRect2（）和BoxPoints（）分别给予碱基编号添加了偏移; 对于所显示的结果，偏移减去回退的计算，使得结果可被容易地比较之后。 理想地，所有各列中的数字（除了牛，OY列）将是相同的; 但作为牛，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


Answer 2:

我碰到这个线程来同时寻找适合的面积最小的边界矩形一个Python的解决方案。

下面是我的实现，这已被证实与Matlab的：
https://stackoverflow.com/a/14675742/1755401



文章来源: Result discrepancy between cv.MinAreaRect2 and ArcGIS (GIS software) . Possible bug?
标签： python debugging opencv geometry
举报

收藏的人(0)

Ta的文章 更多文章
登录 后发表评论
0条评论
还没有人评论过~