Initial commit of Python implementation.

README.txt

@@ -0,0 +1,21 @@
+#
+# Test program to find minimum-area bounding rectangle of a set of 2D points
+# 
+# Currently implemented in Python. I will upload the Matlab version later.
+# 
+#
+# Copyright (c) 2013, David Butterworth, University of Queensland
+# All rights reserved.
+# 
+
+Installation:
+Download all source code from the /python directory to your local computer.
+
+Test the program:
+$ python ./bbox_test.py
+
+The test program includes definitions for some simple polygons.
+You can define your own Nx2 numpy array with your own data.
+
+Tested with Python 2.6.5 on Ubuntu 10.04.4
+Results verified using Matlab

python/bbox_test.py

@@ -0,0 +1,79 @@
+#!/usr/bin/python
+
+# Test program to find minimum-area bounding rectangle
+# 
+# Un-comment one of the arrays to test with a rectangle or 5-point polygon
+# Some solutions to this problem will fail with simple shapes!
+#
+# Tested with Python 2.6.5 on Ubuntu 10.04.4
+# Results verified using Matlab
+
+# Copyright (c) 2013, David Butterworth, University of Queensland
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+#     * Redistributions of source code must retain the above copyright
+#       notice, this list of conditions and the following disclaimer.
+#     * Redistributions in binary form must reproduce the above copyright
+#       notice, this list of conditions and the following disclaimer in the
+#       documentation and/or other materials provided with the distribution.
+#     * Neither the name of the Willow Garage, Inc. nor the names of its
+#       contributors may be used to endorse or promote products derived from
+#       this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+
+from numpy import *
+
+from qhull_2d import *
+from min_bounding_rect import *
+
+if __name__ == "__main__":
+    #
+    # Un-comment one of these shapes below:
+    #
+
+    # Square
+    #xy_points = 10*array([(x,y) for x in arange(10) for y in arange(10)])
+
+    # Random points
+    #xy_points = 100*random.random((32,2))
+    
+    # A rectangle
+    #xy_points = array([ [0,0], [1,0], [1,2], [0,2], [0,0] ])
+
+    # A rectangle, with 5th outlier
+    xy_points = array([ [0,0], [1,0], [1.5,1], [1,2], [0,2], [0,0] ])
+
+    #--------------------------------------------------------------------------#
+
+    # Find convex hull
+    hull_points = qhull2D(xy_points)
+
+    # Reverse order of points, to match output from other qhull implementations
+    hull_points = hull_points[::-1]
+
+    print 'Convex hull points: \n', hull_points, "\n"
+
+    # Find minimum area bounding rectangle
+    (rot_angle, area, width, height, center_point, corner_points) = minBoundingRect(hull_points)
+
+    print "Minimum area bounding box:"
+    print "Rotation angle:", rot_angle, "rad  (", rot_angle*(180/math.pi), "deg )"
+    print "Width:", width, " Height:", height, "  Area:", area
+    print "Center point: \n", center_point # numpy array
+    print "Corner points: \n", corner_points, "\n"  # numpy array
+
+

python/min_bounding_rect.py

@@ -0,0 +1,139 @@
+#!/usr/bin/python
+
+# Find the minimum-area bounding box of a set of 2D points
+#
+# The input is a 2D convex hull, in an Nx2 numpy array of x-y co-ordinates. 
+# The first and last points points must be the same, making a closed polygon.
+# This program finds the rotation angles of each edge of the convex polygon,
+# then tests the area of a bounding box aligned with the unique angles in
+# 90 degrees of the 1st Quadrant.
+# Returns the 
+#
+# Tested with Python 2.6.5 on Ubuntu 10.04.4
+# Results verified using Matlab
+
+# Copyright (c) 2013, David Butterworth, University of Queensland
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+#     * Redistributions of source code must retain the above copyright
+#       notice, this list of conditions and the following disclaimer.
+#     * Redistributions in binary form must reproduce the above copyright
+#       notice, this list of conditions and the following disclaimer in the
+#       documentation and/or other materials provided with the distribution.
+#     * Neither the name of the Willow Garage, Inc. nor the names of its
+#       contributors may be used to endorse or promote products derived from
+#       this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+
+from numpy import *
+import sys  # maxint
+
+def minBoundingRect(hull_points_2d):
+    #print "Input convex hull points: "
+    #print hull_points_2d
+
+    # Compute edges (x2-x1,y2-y1)
+    edges = zeros( (len(hull_points_2d)-1,2) ) # empty 2 column array
+    for i in range( len(edges) ):
+        edge_x = hull_points_2d[i+1,0] - hull_points_2d[i,0]
+        edge_y = hull_points_2d[i+1,1] - hull_points_2d[i,1]
+        edges[i] = [edge_x,edge_y]
+    #print "Edges: \n", edges
+
+    # Calculate edge angles   atan2(y/x)
+    edge_angles = zeros( (len(edges)) ) # empty 1 column array
+    for i in range( len(edge_angles) ):
+        edge_angles[i] = math.atan2( edges[i,1], edges[i,0] )
+    #print "Edge angles: \n", edge_angles
+
+    # Check for angles in 1st quadrant
+    for i in range( len(edge_angles) ):
+        edge_angles[i] = abs( edge_angles[i] % (math.pi/2) ) # want strictly positive answers
+    #print "Edge angles in 1st Quadrant: \n", edge_angles
+
+    # Remove duplicate angles
+    edge_angles = unique(edge_angles)
+    #print "Unique edge angles: \n", edge_angles
+
+    # Test each angle to find bounding box with smallest area
+    min_bbox = (0, sys.maxint, 0, 0, 0, 0, 0, 0) # rot_angle, area, width, height, min_x, max_x, min_y, max_y
+    print "Testing", len(edge_angles), "possible rotations for bounding box... \n"
+    for i in range( len(edge_angles) ):
+
+        # Create rotation matrix to shift points to baseline
+        # R = [ cos(theta)      , cos(theta-PI/2)
+        #       cos(theta+PI/2) , cos(theta)     ]
+        R = array([ [ math.cos(edge_angles[i]), math.cos(edge_angles[i]-(math.pi/2)) ], [ math.cos(edge_angles[i]+(math.pi/2)), math.cos(edge_angles[i]) ] ])
+        #print "Rotation matrix for ", edge_angles[i], " is \n", R
+
+        # Apply this rotation to convex hull points
+        rot_points = dot(R, transpose(hull_points_2d) ) # 2x2 * 2xn
+        #print "Rotated hull points are \n", rot_points
+
+        # Find min/max x,y points
+        min_x = nanmin(rot_points[0], axis=0)
+        max_x = nanmax(rot_points[0], axis=0)
+        min_y = nanmin(rot_points[1], axis=0)
+        max_y = nanmax(rot_points[1], axis=0)
+        #print "Min x:", min_x, " Max x: ", max_x, "   Min y:", min_y, " Max y: ", max_y
+
+        # Calculate height/width/area of this bounding rectangle
+        width = max_x - min_x
+        height = max_y - min_y
+        area = width*height
+        #print "Potential bounding box ", i, ":  width: ", width, " height: ", height, "  area: ", area 
+
+        # Store the smallest rect found first (a simple convex hull might have 2 answers with same area)
+        if (area < min_bbox[1]):
+            min_bbox = ( edge_angles[i], area, width, height, min_x, max_x, min_y, max_y )
+        # Bypass, return the last found rect
+        #min_bbox = ( edge_angles[i], area, width, height, min_x, max_x, min_y, max_y )
+
+    # Re-create rotation matrix for smallest rect
+    angle = min_bbox[0]   
+    R = array([ [ math.cos(angle), math.cos(angle-(math.pi/2)) ], [ math.cos(angle+(math.pi/2)), math.cos(angle) ] ])
+    #print "Projection matrix: \n", R
+
+    # Project convex hull points onto rotated frame
+    proj_points = dot(R, transpose(hull_points_2d) ) # 2x2 * 2xn
+    #print "Project hull points are \n", proj_points
+
+    # min/max x,y points are against baseline
+    min_x = min_bbox[4]
+    max_x = min_bbox[5]
+    min_y = min_bbox[6]
+    max_y = min_bbox[7]
+    #print "Min x:", min_x, " Max x: ", max_x, "   Min y:", min_y, " Max y: ", max_y
+
+    # Calculate center point and project onto rotated frame
+    center_x = (min_x + max_x)/2
+    center_y = (min_y + max_y)/2
+    center_point = dot( [ center_x, center_y ], R )
+    #print "Bounding box center point: \n", center_point
+
+    # Calculate corner points and project onto rotated frame
+    corner_points = zeros( (4,2) ) # empty 2 column array
+    corner_points[0] = dot( [ max_x, min_y ], R )
+    corner_points[1] = dot( [ min_x, min_y ], R )
+    corner_points[2] = dot( [ min_x, max_y ], R )
+    corner_points[3] = dot( [ max_x, max_y ], R )
+    #print "Bounding box corner points: \n", corner_points
+
+    #print "Angle of rotation: ", angle, "rad  ", angle * (180/math.pi), "deg"
+
+    return (angle, min_bbox[1], min_bbox[2], min_bbox[3], center_point, corner_points) # rot_angle, area, width, height, center_point, corner_points
+

python/qhull_2d.py

@@ -0,0 +1,51 @@
+#!/usr/bin/python
+
+# Compute the convex hull of a set of 2D points
+# A Python implementation of the qhull algorithm
+# 
+# Tested with Python 2.6.5 on Ubuntu 10.04.4
+
+# Copyright (c) 2008 Dave (www.literateprograms.org)
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy 
+# of this software and associated documentation files (the "Software"), to deal 
+# in the Software without restriction, including without limitation the rights 
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 
+# copies of the Software, and to permit persons to whom the Software is 
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in 
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS 
+# IN THE SOFTWARE.
+
+from __future__ import division
+from numpy import *
+
+link = lambda a,b: concatenate((a,b[1:]))
+edge = lambda a,b: concatenate(([a],[b]))
+
+def qhull2D(sample):
+    def dome(sample,base): 
+        h, t = base
+        dists = dot(sample-h, dot(((0,-1),(1,0)),(t-h)))
+        outer = repeat(sample, dists>0, 0)
+        if len(outer):
+            pivot = sample[argmax(dists)]
+            return link(dome(outer, edge(h, pivot)),
+                    dome(outer, edge(pivot, t)))
+        else:
+            return base
+    if len(sample) > 2:
+    	axis = sample[:,0]
+    	base = take(sample, [argmin(axis), argmax(axis)], 0)
+        return link(dome(sample, base), dome(sample, base[::-1]))
+    else:
+	return sample
+