The algorithm below fits an ellipse into the polygon (or set of points) in such a way that the area of the ellipse fitted is equal to the area of the polygon.
from __future__ import division
from scipy import *
import pylab as p
def ellipse(ra,rb,ang,x0,y0,Nb=50):
'''ra - major axis length
rb - minor axis length
ang - angle
x0,y0 - position of centre of ellipse
Nb - No. of points that make an ellipse
based on matlab code ellipse.m written by D.G. Long,
Brigham Young University, based on the
CIRCLES.m original
written by Peter Blattner, Institute of Microtechnology,
University of
Neuchatel, Switzerland, blattner@imt.unine.ch
'''
xpos,ypos=x0,y0
radm,radn=ra,rb
an=ang
co,si=cos(an),sin(an)
the=linspace(0,2*pi,Nb)
X=radm*cos(the)*co-si*radn*sin(the)+xpos
Y=radm*cos(the)*si+co*radn*sin(the)+ypos
return X,Y
def _makexy(scale):
x=reshape(array(range(scale*2)*(scale*2)),(scale*2,scale*2))
return x,x.T
def _inside2(x,y,px,py):
'''Check if point (x,y) inside polygon deifined by (px,py).
From: http://www.ariel.com.au/a/python-point-int-poly.html
'''
poly=zip(px,py)
n = len(poly)
inside =False
p1x,p1y = poly[0]
for i in range(n+1):
p2x,p2y = poly[i % n]
if y > min(p1y,p2y):
if y <= max(p1y,p2y):
if x <= max(p1x,p2x):
if p1y != p2y:
xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
if p1x == p2x or x <= xinters:
inside = not inside
p1x,p1y = p2x,p2y
return inside
def _polyMask(X,Y,rs):
'''Make binary mask of polygon. Values of one
in 2D array pd represents points enclosed by the polygon.'''
pd=zeros((rs*2,rs*2))
for i in range(pd.shape[0]):
for j in range(pd.shape[1]):
if _inside2(i,j,X,Y)==True:
pd[i,j]=1.0
return pd
def ellfit(X,Y,rs=50,showFig=True):
'''Fit and ellipse to a polygon defined by X,Y in such a way
that area of the ellispe fitted is the same as area of the polygon.
Based on ellfit.pro by D.G. Long, Brigham Young University.
Conditions:
1. The polygon must not be self intersecting.
2. Polygon must be closed (last point must match the first one).
3. Polygon must be centered in (0,0)
4. Polygon must be within rectangle <-1,1>x<-1,1>
'''
#Xmax,Ymax=abs(X).max(),abs(Y).max()
#X=X/Xmax
#Y=Y/Ymax
pi2=2*pi
rs2=2*rs
if showFig:
fig = p.figure(figsize=(5,5))
ax = fig.add_subplot(111)
X=X*rs+rs
Y=Y*rs+rs
xa,ya=_makexy(rs)
pd=_polyMask(X,Y,rs)
if showFig:
p.plot(X,Y,"rs-",ms=15)
for i in range(pd.shape[0]):
for j in range(pd.shape[1]):
if pd[i,j]==1.0:
pass
p.plot([i],[j],'g.',ms=10)
m00=pd.sum()
m10=(xa*pd).sum()
m01=(ya*pd).sum()
xm,ym=m10/m00, m01/m00
u20=(xa*xa*pd).sum()-xm*m10
u02=(ya*ya*pd).sum()-ym*m01
u11=(xa*ya*pd).sum()-ym*m10
tn2=2*u11/(u20-u02)
ang=arctan(tn2)/2.0
if u02>u20:
ang=ang-pi/2
xmean,ymean=X.mean(),Y.mean()
Xr = cos(ang) * (X-xmean) - sin(ang) * (Y-ymean) +xmean
Yr = sin(ang) * (X-xmean) + cos(ang) * (Y-ymean) +ymean
t=_polyMask(Xr,Yr,rs)
m00=t.sum()
m10=(xa*t).sum()
m01=(ya*t).sum()
xm,ym=m10/m00, m01/m00
u20=(xa*xa*t).sum()-xm*m10
u02=(ya*ya*t).sum()-ym*m01
u11=(xa*ya*t).sum()-ym*m10
aa=((16*u20**3)/(pi2*u02))**(1/8)
bb=((16*u02**3)/(pi2*u20))**(1/8)
##ellispe params
a,b=max([aa,bb]),min([aa,bb])
b=t.sum()/(pi*a)
ecc=sqrt(a**2-b**2)
theta=ang
ellArea=pi*a*b
print "Ploygon_area - ellispe_area = ",round(ellArea-t.sum(),3)
#print time.time()-t1
if showFig:
p.grid(True)
p.show()
#xee,yee=b/rs*cos(0),a/rs*sin(0)
#b3,a3=xee*max([Xmax,Ymax])/cos(theta),yee*min([Xmax,Ymax])/sin(theta)
#print b3/a3, b/a
return a/rs,b/rs,ecc,theta
def elltest(scale=0.8,off=0.2):
#generate an example, random, non-self-intersecting polygon.
#This is done by first generating
#it in polar coordinates and than translating it
#to cartesian.
Theta1,R1=linspace(0,2*pi,30),rand(30)*scale+off
X1,Y1=R1*cos(Theta1),R1*sin(Theta1)
X1=append(X1,X1[0])
Y1=append(Y1,Y1[0])
p.plot(X1,Y1,".-",ms=10)
a2,b2,ecc2,alpha2=ellfit(X1,Y1,showFig=False)
Xe,Ye=ellipse(b2,a2,-alpha2,X1.mean(),Y1.mean(),Nb=40)
p.plot(Xe,Ye,"r.-")
p.grid(True)
p.show()
pass
if __name__ == '__main__':
elltest()
Examples
Note 1
The algorithm has the following limitations:
- The polygon must not be self intersecting.
- Polygon must be closed (last point must match the first one).
- Polygon must be centered in (0,0)
- Polygon must be within rectangle <-1,1>x<-1,1>
Note 2
This algorithm could be modified to overcome the above (Note 1) limitations. However, at the very moment I don't have time to do this. Hope this is not a problem.
Note 3
Searching the Internet for an ellipse fitting script for Python I found
fitellipse.py. It is very useful, however, sometimes I got an Exception that its not possible to fit an ellipse. The algorithm presented here seems to work good. In case the ellipse cannot be fitted try to increase rs attribute of ellfit function. The default value is 50 (rs=50).
Note 3
This script is in
ellipsefit.py.
