Selecting a maximal area quadrilateral from four random points

To calculate the area of a quadrilateral with the points coordinates, the points are supposed to be taken clockwise or counter clock wise.When the shoelace formula is applied on each of  4!=24 tuple of four points, different values can be obtained:

For example, the area of the square (on top left) with the points taken clock wise, is equal to 4. The 24 ways to calculate the area yield only two values 0 or 4. The quadrilateral of area equal to 0 is plotted on top right (red), this quadrilateral is self intersecting. The trapeze area (down left) is equal to 3, the 24 ways to calculate the area yied three values :0, 1, 3 and the maximal value of the area corresponds to a non intersecting quadrilateral (blue, down right).

Quadrilaterals were build from randomly chosen points, their area was calculated by taking the original points order:

Some of them are self intersecting (first, top left) or degenerated (the 10th). By searching the order of the points giving the maximal area, non self-intersecting quadrilaterals can be found:

If both initial and final quadrilaterals are overlayed:

Compare from the previous post, the function to calculate the area was modified http://people.virginia.edu/~ll2bf/docs/various/polyarea.html).

Play with the code on cloud.sagemath or download code.

Update (2014-04):

A quadrilateral of maximal area is found using brute force (by looking for the permutations of four points). Since cyclic or circular permutations of four points yield quadrilaterals with the same area, the search of maximal area should be performed not on a set of permutations of points but on a set of unique cyclic permutations of four points (a "necklace") possibly with Sawada's algorithm

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

view raw getquadrilateral.py hosted with ❤ by GitHub