Skip to content

Instantly share code, notes, and snippets.

@jeanpat
Created May 27, 2013 10:13
Show Gist options
  • Save jeanpat/5656335 to your computer and use it in GitHub Desktop.
Save jeanpat/5656335 to your computer and use it in GitHub Desktop.
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 31 10:40:59 2013
@author: Jean-Patrick Pommier
http://dip4fish.blogspot.fr/2013/01/selecting-simple-quadrilateral-from.html
simpoly :
http://people.virginia.edu/~ll2bf/docs/various/polyarea.html
"""
import random, math
import itertools as it
from itertools import cycle
from collections import defaultdict
import numpy as np
from scipy.misc import comb
from matplotlib import cm
from matplotlib import pyplot as plt
from matplotlib.patches import Polygon
def simpoly(x,y):
"""
A function that calculates the area of a 2-D simple polygon (no matter concave or convex)
Must name the vertices in sequence (i.e., clockwise or counterclockwise)
Square root input arguments are not supported
Formula used: http://en.wikipedia.org/wiki/Polygon#Area_and_centroid
Definition of "simply polygon": http://en.wikipedia.org/wiki/Simple_polygon
Input x: x-axis coordinates of vertex array
y: y-axis coordinates of vertex array
Output: polygon area
"""
ind_arr = np.arange(len(x))-1 # for indexing convenience
s = 0
for ii in ind_arr:
s = s + (x[ii]*y[ii+1] - x[ii+1]*y[ii])
return abs(s)*0.5
#def quadAreaShoelace(A,B,C,D):
# x1 = A[0]
# y1 = A[1]
# x2 = B[0]
# y2 = B[1]
# x3 = C[0]
# y3 = C[1]
# x4 = D[0]
# y4 = D[1]
# #sign +/- if direct/indirect quadrilateral
#
# return 0.5*abs(x1*y2+x2*y3+x3*y1-x2*y1-x3*y2-x4*y3-x1*y4)
def maxAreaQuad(A, B, C, D):
allquads = it.permutations((A, B, C, D))
quads = [q for q in allquads]
#print len(quads)
dicQuads = defaultdict(list)
for quad in quads:
#print 'quad in sholace:', quad
#area = quadAreaShoelace(*quad)
x = [p[0] for p in quad]
y = [p[1] for p in quad]
area = simpoly(x,y)
dicQuads[area].append(quad)
allkeys = dicQuads.keys()
#print type(allkeys), allkeys
#sortedkeys = sorted[allkeys]
allkeys.sort()
areamax = allkeys[-1]
return areamax,allkeys,dicQuads[areamax][0]#sorted(set(allkeys)),areamax,
def plotQuad(quadrilateral,col='g',alph=0.01):
p = Polygon( quadrilateral, alpha=alph, color=col )
plt.gca().add_artist(p)
if __name__ == "__main__":
Nb_points = 5
somepoints = [(1, 0), (0, 0), (2, 2), (1, 3), (0, 2),(1,1),(20,20)]
A = 0,0
B = 2,0
C = 1,2
D = 0,2
E = 2,2
F = 1,1
square = A, B, E, D
selfsq = A, B, D, E
trapez = A, B, C, D
selftz = A, C, B, D
geofig = square, selfsq,trapez, selftz
n=1
fig1 = plt.figure()
for f in geofig:
# print f
# print maxAreaQuad(*f)
ax=plt.subplot(2,2,n,frameon = True,xticks = [0,1,2], yticks = [0,1,2])#
for pt in f:
ax.scatter(*pt,c='blue',s=50)
n =n +1
x = [p[0] for p in f]
y = [p[1] for p in f]
area1 = simpoly(x,y)
plotQuad(f,alph=0.5, col='r')
areamax,allareas,cfig = maxAreaQuad(*f)
plotQuad(cfig,alph=0.1, col='b')
plt.title(str(area1)+' : '+str(allareas), fontsize=8)
fig2 = plt.figure()
random.seed(a=52521)
for n in range(1,26):
points=[(random.randint(0,20),random.randint(0,20)) for i in range(0,4)]
#print 'points ', len(points)
quads=[i for i in it.combinations(points,4)]
#print len (quads)
for q in quads:
#print q
#print maxAreaQuad(*q)
ax1=plt.subplot(5,5,n,frameon = True,xticks = [], yticks = [])#
for pt in q:
ax1.scatter(*pt,c='blue',s=50)
x = [p[0] for p in q]
y = [p[1] for p in q]
area1 = simpoly(x,y)
areamax,allareas,maxquad = maxAreaQuad(*q)
# inputquad_area = quadAreaShoelace(*q)
# outputquad_area = quadAreaShoelace(*maxquad)
plotQuad(q,alph=1, col='r')
plotQuad(maxquad,alph=0.5, col='b')
#shoelace implementation
# sh_area_in= str(int(inputquad_area))
# sh_area_out = str(int(outputquad_area))
# sh_area_max = str(int(areamax))
#simple area
s_area_in = simpoly(x,y)
xx=[p[0] for p in maxquad]
yy=[p[1] for p in maxquad]
s_area_out = simpoly(xx,yy)
#ti = sh_area_in+'('+sh_area_out+'='+sh_area_max+')'
ti = '['+str(s_area_in)+']'#+ti
ti = ti+'['+str(s_area_out)+']'
#plt.title(ti, fontsize=8)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment