Result discrepancy between cv.MinAreaRect2 and Arc

Question

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>

Answer 1

I came across this thread while looking for a Python solution for a minimum-area bounding rectangle.

Here's my implementation, which has been verified with Matlab:
https://stackoverflow.com/a/14675742/1755401

Answer 2

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