getquadrilateral.py

# -*- 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()