Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 27 additions & 17 deletions python/bbox_test.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#!/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!
#
Expand Down Expand Up @@ -46,34 +46,44 @@
#

# Square
#xy_points = 10*array([(x,y) for x in arange(10) for y in arange(10)])
# 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))
# xy_points = 100*random.random((32,2))

# A rectangle
#xy_points = array([ [0,0], [1,0], [1,2], [0,2], [0,0] ])
# 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] ])
v = [
[282.15683092396523, 223.11512640953208],
[282.6740563777786, 200.61399296154428],
[248.98051148730354, 199.7099293962484],
[248.46328600395807, 217.4898437932211],
[229.3259436222214, 273.2403900872115],
[261.3200334037592, 284.18959062917554],
[282.15683092396523, 223.11512640953208],
]

xy_points = array(v)

#--------------------------------------------------------------------------#
# --------------------------------------------------------------------------#

# 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"
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


(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
135 changes: 82 additions & 53 deletions python/min_bounding_rect.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@

# 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 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
# Returns the
#
# Tested with Python 2.6.5 on Ubuntu 10.04.4
# Results verified using Matlab
Expand Down Expand Up @@ -40,100 +40,129 @@
# POSSIBILITY OF SUCH DAMAGE.

from numpy import *
import sys # maxint
import sys # maxsize
import math


def minBoundingRect(hull_points_2d):
#print "Input convex hull points: "
#print 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
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
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
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
# 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) ):
min_bbox = (
0,
sys.maxsize,
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
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
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
# 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
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 )
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 )
# 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
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
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
# 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
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

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
48 changes: 24 additions & 24 deletions python/qhull_2d.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,50 +2,50 @@

# 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
# 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
# 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
# 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]))
link = lambda a, b: concatenate((a, b[1:]))
edge = lambda a, b: concatenate(([a], [b]))


def qhull2D(sample):
def dome(sample,base):
def dome(sample, base):
h, t = base
dists = dot(sample-h, dot(((0,-1),(1,0)),(t-h)))
outer = repeat(sample, dists>0, 0)
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)))
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)
axis = sample[:, 0]
base = take(sample, [argmin(axis), argmax(axis)], 0)
return link(dome(sample, base), dome(sample, base[::-1]))
else:
return sample

return sample