Несоответствие результатов между cv.MinAreaRect2 и ArcGIS (программное обеспечение ГИС). Возможная ошибка?

У меня есть набор точек, полученных из многоугольника. Я тестирую несколько решений, чтобы получить минимальную площадь или прямоугольник. В качестве эталона я использую 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 находит ограниченный прямоугольник минимальной области для двумерной точки, созданной построением выпуклой оболочки для множества и применением техники вращающихся суппортов к корпусу.

 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 в шейп- box_vtx для просмотра в ArcGIS и сравниваю его с минимальной границей геометрии (Data Management), я вижу эту разницу, как показано на следующем рисунке, где:

  • красный = граница многоугольника
  • синий = прямоугольник с минимальной площадью от 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> 

Я не смотрел код OpenCV, но кажется вероятным, что большие x, y продукты вычитаются из других больших x, y продуктов в его вычислениях. Смещения в ваших значениях x используют около 19 бит, а значения у около 23, поэтому такое вычитание может привести к потере около 42 бит из 53, которые переносят типичные числа двойной точности. (См. Плавающие точки в википедии.) Размер и форма желтого и черного прямоугольников выглядят разумно, но отображаемая ширина и длина (10.285 …, 18.335 …) на 1% отличаются от (10.393, 18.037), ширины и длины, которые возникли в соответствующем вопросе .

Короче говоря, 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). Значение выражения карты – это список переведенных точек. Обратите внимание, что может быть предпочтительнее установить cx = int(sum(zip(*s)[0])/n) и cy = int(sum(zip(*s)[1])/n) так что точки решетки остаются точки решетки, если это то, что вы вычисляете. Благоприятный эффект удаления больших смещений остается, и во время перевода будет происходить меньшее округление.

Примечание. Я провел тесты с смещениями, вычитаемыми из набора данных и из корпуса (предположим, что корпус показан в вопросе № 13542855 и, возможно, в вопросе № 13553884 ) и получил несогласованные результаты, которые не согласуются с результатами метода, который я дал в № 13542855. Таким образом, трудно определить, где возникает проблема, исходя из ваших чисел. Ниже показан более простой тестовый пример. Этот тест дает понять, что большие смещения вызывают низкую точность. Возможно, вы можете запускать аналогичные тестовые данные через ArcGIS. Ниже приведены результаты теста. MinAreaRect2 () и BoxPoints () получили базовые номера с добавленными смещениями; для отображаемых результатов смещение вычитается после вычисления, поэтому результаты можно легко сравнивать. В идеале все числа в каждом столбце (кроме столбцов ox, oy) будут одинаковыми; но по мере роста ox, они начинают различаться и вскоре становятся почти случайными.

  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 

Я наткнулся на эту тему, ища решение Python для прямоугольника с ограниченной границей.

Вот моя реализация, которая была проверена с помощью Matlab:
https://stackoverflow.com/a/14675742/1755401