cv.MinAreaRect2和ArcGIS（GIS软件）之间的结果差异。可能的错误？ [英] Result discrepancy between cv.MinAreaRect2 and ArcGIS (GIS software) . Possible bug?
问题描述
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， （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.4 495758876,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（）。
函数cv.MinAreaRect2通过为集合构建凸包来为2D点集合找到
最小区域的外接矩形，并将
应用于船体的旋转卡尺技术。
import cv
＃（x，y） - 框的中心点
＃（w， h） - 框的宽度和高度
＃theta - 旋转角
（（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中查看并与：
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）
＃获取多边形边界的点（上面的点）
points = []
for x in xrange（pts.GetPointCount（））：
points.append（（pts.GetX（p），pts.GetY（p）））
＃Convex Hull
CH1 = geometry.ConvexHull
＃我没有找到一个方法来限制点
print CH1（）
＃使用openCV
cvxHull = cv.ConvexHull2（points， cv.CreateMemStorage（），return_points = True）
print cvxHull
< cv2.cv.cvseq对象在0x00000000067CCF90>
我没有看过OpenCV代码，在其计算中似乎可能会从其他大型x，y产品中减去大量的x，y产品。 x值中的偏移量约为19位，y的偏移约为23，因此这种减法可能导致典型双精度数字携带的53位丢失约42位。 （请参阅维基百科中的浮点）。黄色和黑色的大小和形状矩形看起来很合理但显示的宽度和长度（10.285 ...，18.335 ...）与（10.393，18.037）1％不同，宽度和长度出现在相关问题。
简而言之，OpenCV可能在 BoxPoints（）。
要检查或测试的事情开发包括：
•显示和打印OpenCV的凸包计算点
•在图形上，显示OpenCV和ArcGIS矩形的中心，并打印距离这些中心的距离显示的矩形的角落
•使用已翻译的数据集（见下文）重新运行计算和图形，以减少下溢ef相互减去大数字的方式
在OpenCV计算之前翻译数据集可以将丢失的位数减少一半。通过以下代码生成一组新的点，并使用新的数据集进行计算：
s = points＃积分=原始数据集
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）
＃现在点=翻译数据集
在上述代码中， zip（* s）解压缩（x ，y）指向两个列表，x值列在zip（* s）[0]中，y值列在zip（* s）[1]中。因此（cx，cy）表示s中列出的点的质心。地图表达式将函数应用于可迭代的元素。该函数返回由（-cx，-cy）翻译的点。地图表达式的值是一个翻译点列表。注意，最好设置 cx = int（sum（zip（* s）[0]）/ n）和 cy = int sum（zip（* s）[1]）/ n），使格点保持格点，如果这是你计算的。删除大偏移量的有益效果仍然存在，而在翻译过程中会发生较少的舍入。
注意 - 我运行了从数据中减去的偏移量的测试设置和从船体（假设船体如问题＃13542855 ，也许在问题＃13553884 ），并得到不一致的结果与我在＃13542855中给出的方法的结果不一致。因此，根据您的数字难以确定出现问题的位置。一个更简单的测试用例如下所示。该测试清楚地表明，大的偏移导致精度差。也许您可以通过ArcGIS运行类似的测试数据。以下是测试结果的一半。 MinAreaRect2（）和BoxPoints（）被赋予了带有偏移量的基数;对于显示结果，计算后减去偏移量，以便可以轻松比较结果。理想情况下，每列（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}'。格式（
ox，oy，x-ox，y-oy，w ，h，th））
ox，oy = 4 ** i，5 ** i
print'\\\
ox oy T x T y T x T y T x T y在框中的（ox，oy，box）框中的T x T y'
print'{：10d} {：10d}'。format（ox，oy）
for p in box：
print'{：9.4f} {：9.4f}'。format（p [0] -ox，p [1] -oy），
print
I have a set of points derived from a polygon. I am testing several solution in order to obtain the minimum area or the rectangle. As benchmark I am using 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)]
One solution is use cv.MinAreaRect2() of OpenCV.
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)
when I convert box_vtx in a shapefile in order to view in ArcGIS and compare with Minimum Bounding Geometry (Data Management)I can see this difference as shown in the following picture,
where:
- red = border of polygon
- blue = rectangle with minimum area from ArGIS (10.1)
- yellow and black = rectangle with minimum area from OpenCV
Working with OpenCV compared with solution proposed in this 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>
I haven't looked at OpenCV code, but it seems likely that large x,y products are being subtracted from other large x,y products in its calculations. The offsets in your x values use about 19 bits, and those of y about 23, so such subtraction can result in a loss of about 42 bits from the 53 that typical double precision numbers carry. (See floating point in wikipedia.) The size and shape of the yellow and black rectangle look reasonable but the displayed width and length, (10.285..., 18.335...) are 1% different from (10.393, 18.037), the width and length that appeared in a related question.
In short, OpenCV may have a rounding, underflow, or overflow problem in MinAreaRect2() or in BoxPoints().
Things to check or tests to develop include:
• Display and print the points of OpenCV's convex hull computation
• On the graph, display the centers of OpenCV's and ArcGIS's rectangles, and print distances from those centers to the corners of the displayed rectangles
• Rerun the computations and graph with a translated data set (see below), to reduce the underflow effect of subtracting large numbers from each other
Translating the data set before the OpenCV computation can cut the number of bits lost in half. Produce a new set of points via code like the following, and try the computation with the new data set:
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
In the above code, zip(*s) unzips the list of (x,y) points into two lists, with the x values listed in zip(*s)[0] and the y values listed in zip(*s)[1]. Thus (cx,cy) represents the center of mass of points listed in s. The map expression applies a function to elements of an iterable. The function returns a point translated by (-cx,-cy). The value of the map expression is a list of translated points. Note, it might be preferable to set cx = int(sum(zip(*s)[0])/n) and cy = int(sum(zip(*s)[1])/n) so that lattice points remain lattice points, if that is what you are computing with. The beneficial effect of removing large offsets remains, and less rounding will occur during translation.
Note – I ran tests with offsets subtracted from the data set and from the hull (supposing the hull is as shown in question #13542855 and perhaps in question #13553884) and got inconsistent results that do not agree well with results from the method I gave in #13542855. Thus it is difficult to pin down where the problem is occurring based on your numbers. A simpler test case is shown below. This test makes it clear that large offsets cause poor accuracy. Perhaps you could run similar test data via ArcGIS. Following are half of the results from the test. MinAreaRect2() and BoxPoints() were given base numbers with offsets added; for the displayed results, the offset is subtracted back off after the calculation so that results can be compared easily. Ideally, all of the numbers in each column (except ox,oy columns) would be the same; but as ox,oy increase, they begin to differ and soon become nearly random.
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
Here is the code that produced the data shown above:
#!/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
这篇关于cv.MinAreaRect2和ArcGIS（GIS软件）之间的结果差异。可能的错误？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
