Skip to content

Instantly share code, notes, and snippets.

@kchr
Created August 31, 2015 11:40
Show Gist options
  • Save kchr/77a0ee945e581df7ed25 to your computer and use it in GitHub Desktop.
Save kchr/77a0ee945e581df7ed25 to your computer and use it in GitHub Desktop.
Python implementation: Convex hull + Minimal bounding rectangle

What is this?

Python proof-of-concept implementation of two geomapping algorithms

Requirements

This code use the following packages:

  • numpy
  • matplotlib (optional, only for creating graphs)

How to use

$ python ./test.py

#!/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
import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection
import matplotlib as mpl
def poly_plot(verts):
# Generate data. In this case, we'll make a bunch of center-points and generate
# verticies by subtracting random offsets from those center-points
numpoly, numverts = 100, 4
centers = 100 * (np.random.random((numpoly,2)) - 0.5)
offsets = 10 * (np.random.random((numverts,numpoly,2)) - 0.5)
#verts = centers + offsets
#verts = np.swapaxes(verts, 0, 1)
# In your case, "verts" might be something like:
# verts = zip(zip(lon1, lat1), zip(lon2, lat2), ...)
# If "data" in your case is a numpy array, there are cleaner ways to reorder
# things to suit.
# Color scalar...
# If you have rgb values in your "colorval" array, you could just pass them
# in as "facecolors=colorval" when you create the PolyCollection
z = np.random.random(numpoly) * 500
fig, ax = plt.subplots()
# Make the collection and add it to the plot.
coll = PolyCollection(verts, array=z, cmap=mpl.cm.jet, edgecolors='none')
ax.add_collection(coll)
ax.autoscale_view()
# Add a colorbar for the PolyCollection
fig.colorbar(coll, ax=ax)
plt.show()
#!/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
import math
import sys
#from barisal_coords import coords
from borgona_coords import coords
#from bhola_coords import coords
#from bangalore_coords import coords
from numpy import *
from qhull_2d import *
from min_bounding_rect import *
from polyplot import poly_plot
grid = []
for lon, lat in coords:
grid.append([lon, lat])
xy_points = array(grid)
# Find convex hull
hull_points = qhull2D(xy_points)
print len(hull_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)
# Verbose output of return data
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
# Draw a nice graph with the new shape
poly_plot(array([corner_points]))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment