A Python Script to Fit an Ellipse to Noisy Data

Problem statement

Given a set of noisy data which represents noisy samples from the perimeter of an ellipse, estimate the parameters which describe the underlying ellipse.

Discussion

There are two general ways to fit an ellipse: algebraic and geometric approaches. In an algebraic approach, the parameters for an algebraic description of an ellipse are fit subject to constraints which guarantee the parameters result in an ellipse. In the geometric approach,  characteristics of the ellipse are fit.

The code snippet below uses a method described by Yu, Kulkarni & Poor. The location of the foci and the length of the line segments from the foci to a point on the perimeter of the ellipse are found through an optimization problem. Because the fitting objective is not convex and has a minimum at infinity, a penalty cost is added to prevent the foci from wandering off.

Code

'''
Script to fit an ellipse to a set of points.
- The ellipse is represented by the two foci and the length of a 
     line segment which is drawn from the foci to the 
     point where the ellipse intersects the minor axis.
    
- Fitting algorithm from Yu, Kulkarni & Poor

'''

__author__ = 'Ed Tate'
__email__  = 'edtategmail-dot-com'
__website__ = 'exnumerus.blogspot.com'
__license__ = 'Creative Commons Attribute By - http://creativecommons.org/licenses/by/3.0/us/'''

####################################################
# create ellipse with random noise in points
from random import uniform,normalvariate
from math import pi, sin, cos, exp, pi, sqrt
from openopt import NLP
from numpy import *
from numpy import linalg as LA
import matplotlib.pylab as pp

def gen_ellipse_pts(a,foci1,foci2,
                    num_pts=200, angles=None,
                    x_noise = None, y_noise=None):

    '''
       Generate points for an ellipse given
          the foci, and
          the distance to the intersection of the minor axis and ellipse.
    
       Optionally, 
          the number of points can be specified,
          the angles for the points wrt to the centroid of the ellipse, and 
          a noise offset for each point in the x and y axis.
    '''
    c = (1/2.0)*LA.norm(foci1-foci2)
    b = sqrt(a**2-c**2)
    x1 = foci1[0]
    y1 = foci1[1]
    x2 = foci2[0]
    y2 = foci2[1]
    if angles is None:
        t = arange(0,2*pi,2*pi/float(num_pts))
    else:
        t = array(angles)
            
    ellipse_x = (x1+x2)/2 +(x2-x1)/(2*c)*a*cos(t) - (y2-y1)/(2*c)*b*sin(t)
    ellipse_y = (y1+y2)/2 +(y2-y1)/(2*c)*a*cos(t) + (x2-x1)/(2*c)*b*sin(t)
    try:
        # try adding noise to the ellipse points
        ellipse_x = ellipse_x + x_noise
        ellipse_y = ellipse_y + y_noise
    except TypeError:
        pass
    return (ellipse_x,ellipse_y)

####################################################################

# setup the reference ellipse

# define the foci locations
foci1_ref = array([2,-1])
foci2_ref = array([-2,1])
# pick distance from foci to ellipse
a_ref = 2.5

# generate points for reference ellipse without noise
ref_ellipse_x,ref_ellipse_y = gen_ellipse_pts(a_ref,foci1_ref,foci2_ref)

# generate list of noisy samples on the ellipse
num_samples = 1000
angles = [uniform(-pi,pi) for i in range(0,num_samples)]
sigma = 0.2
x_noise = [normalvariate(0,sigma) for t in angles]
y_noise = [normalvariate(0,sigma) for t in angles]
x_list,y_list = gen_ellipse_pts(a_ref,foci1_ref,foci2_ref,
                                angles  = angles,
                                x_noise = x_noise,
                                y_noise = y_noise)

point_list = []
for x,y in zip(x_list,y_list):
    point_list.append(array([x,y]))    

# draw the reference ellipse and the noisy samples    
pp.figure()
pp.plot(x_list,y_list,'.b', alpha=0.5)
pp.plot(ref_ellipse_x,ref_ellipse_y,'g',lw=2)
pp.plot(foci1_ref[0],foci1_ref[1],'o')
pp.plot(foci2_ref[0],foci2_ref[1],'o')

#####################################################

def initialize():
    '''
    Determine the initial value for the optimization problem.
    '''
    # find x mean
    x_mean = array(x_list).mean()
    # find y mean
    y_mean = array(y_list).mean()
    # find point farthest away from mean
    points = array(zip(x_list,y_list))
    center = array([x_mean,y_mean])
    distances = zeros((len(x_list),1))
    for i,point in enumerate(points):
        distances[i,0]=LA.norm(point-center)
    ind = where(distances==distances.max())
    max_pt = points[ind[0],:][0]
    # find point between mean and max point
    foci1 = (max_pt+center)/2.0
    # find point opposite from 
    foci2 = 2*center - max_pt
    return [distances.max(), foci1[0],foci1[1],foci2[0],foci2[1]]


def objective(x):
    '''
    Calculate the objective cost in the optimization problem.
    '''
    foci1 = array([x[1],x[2]])
    foci2 = array([x[3],x[4]])
    a     = x[0]
    n = float(len(point_list))
    _lambda =0.1
    _sigma = sigma
    sum = 0
    for point in point_list:
        sum += ((LA.norm(point-foci1,2)+LA.norm(point-foci2,2)-2*a)**2)/n
    sum += _lambda*ahat_max*_sigma*exp((a/ahat_max)**4)
    return sum

# solve the optimization problem
x0 = initialize()
ahat_max = x0[0]
print x0
p = NLP(objective, x0)
r = p.solve('ralg')
print r.xf

# get the results from the optimization problem
xf = r.xf
# unload the specific values from the result vector
foci1 = array([xf[1],xf[2]])
foci2 = array([xf[3],xf[4]])
a     = xf[0]

# reverse the order of the foci to get closest to ref foci
if LA.norm(foci1-foci1_ref)>LA.norm(foci1-foci2_ref):
    _temp = foci1
    foci1 = foci2
    foci2 = _temp

####################################################
# plot the fitted ellipse foci
pp.plot([foci1[0]],[foci1[1]],'xk')
pp.plot([foci2[0]],[foci2[1]],'xk')

# plot a line between the fitted ellipse foci and the reference foci
pp.plot([foci1[0],foci1_ref[0]],[foci1[1],foci1_ref[1]],'m-')
pp.plot([foci2[0],foci2_ref[0]],[foci2[1],foci2_ref[1]],'m-')

# plot fitted ellipse
(ellipse_x,ellipse_y) = gen_ellipse_pts(a,foci1,foci2,num_pts=1000)  
pp.plot(ellipse_x,ellipse_y,'r-',lw=3,alpha=0.5)

# scale the axes for a square display
x_max = max(x_list)
x_min = min(x_list)
y_max = max(y_list)
y_min = min(y_list)

box_max = max([x_max,y_max])
box_min = min([x_min,y_min])
pp.axis([box_min, box_max, box_min, box_max])

pp.show()

---

References

• Robust Fitting of Ellipses and Spheroids by Yu, Kulkarni, & Poor
• Least-Squares Fitting of Circles and Ellipses
• http://www.cs.cornell.edu/cv/OtherPdf/Ellipse.pdf
• http://www.math.temple.edu/~renault/ellipses.html
• Numerically Stable Direct Least Squares Fitting of Ellipses
• Notes on plotting an ellipse from foci and distance

---

Testing Configuration

• PythonXY 2.5.6.2
  • Matplotlib 1.0.0
• OpenOpt 0.32

This work is licensed under a Creative Commons Attribution By license.

A Python Recipe to Interactively Identify Elements in Matplotlib

Problem Statement

You need to interactively identify a point in a Matplotlib plot.

Solution

This solution requires three parts. First, a function is defined which is executed when a plot element is selected. Second, this function is connected to the figure. Finally, when a point is plotted, a picker radius is added which defines how big of an area should react to the pick event. One way to get a unique response from each point is to add an attribute to each point. In this example, a string is added with information about the point location. More complex responses can be created by considering an attribute like this in ‘onpick’ function.

__author__ = 'Ed Tate'
__email__  = 'edtate<at>gmail-dot-com'
__website__ = 'exnumerus.blogspot.com'
__license__ = 'Creative Commons Attribute By - http://creativecommons.org/licenses/by/3.0/us/'''

import matplotlib.pylab as p
import random

# define what should happen when a point is picked
def onpick(event):
    thisline = event.artist
    name = thisline.name
    print name

# create the figure and add the handler which reacts to a pick event
fig = p.figure()
fig.canvas.mpl_connect('pick_event',onpick)

# create random set of points to plot
x_pts = [random.uniform(0,1) for i in range(0,500)]
y_pts = [random.uniform(0,1) for i in range(0,500)]

for x,y in zip(x_pts,y_pts):
    pt, = p.plot(x,y,'.b',picker=3)
    pt.name = 'Point located at (%f,%f)' % (x,y)
    

p.show()

When this script is executed, picking a point results in a message which shows which point was chosen:

Point located at (0.481866,0.661643)

---

References

• Matplotlib documentation and example
• Matplotlib pick event demo

---

Testing Conditions

• PythonXY 2.5.6.2
• Matplotlib 1.0.0

This work is licensed under a Creative Commons Attribution By license.

Rough Draft: How to generate 2D Delaunay Triangulations using Qhull and Python

Delaunay triangulations partition a space into regions. Delaunay triangulations can be useful for interpolation and visualizing spatial data.

Qhull is a program which can generate tesselations, convex hulls, and vonoroi diagrams from a set of points. This program is available as a precompiled executable and source code. By interfacing to the command line version of this program, a Delaunay triangulation can be generated.

Sample CodE

To use this code, download Qhull and copy the ‘qhull.exe’ to the working directory. Make sure the code has read, write, and execute privileges in the working directory. Make sure there are not files named ‘data.txt’ or ‘results.txt’ which need to be preserved.

'''
This module generates 2D tesselations for a set of points.
The vornoi cells can be filtered by supplying a value associated with 
   each node.
'''

__author__ = 'Ed Tate'
__email__  = 'edtate-at-gmail-dot-com'

def delaunay2D(xpt,ypt,cpt=None,threshold=0):
    
    if cpt is None:
        cpt = [-1 for x in xpt]    
        
    # write the data file
    pts_filename = 'data.txt'
    pts_F = open(pts_filename,'w')
    #print pts_F
    pts_F.write('2 # this is a 2-D input set\n')
    pts_F.write('%i # number of points\n' % len(xpt))
    for i,(x,y) in enumerate(zip(xpt,ypt)):
        pts_F.write('%f %f # data point %i\n' % (x,y,i))
    pts_F.close()

    # trigger the shell command
    import subprocess

    p = subprocess.Popen('qhull TI data.txt TO results.txt d i Qc Qt Qbb', shell=True)
    p.wait()

    # open the results file and parse results
    results = open('results.txt','r')
    print results

    # get 'i' results
    data = results.readline()
    tri_list = []
    numLines = int(data)
    print numLines
    for i in range(numLines):
        # load each triplet of indexes to x,y points
        data = results.readline()
        idx1,idx2,idx3,dummy = data.split(' ')
        idx1 = int(idx1)
        idx2 = int(idx2)
        idx3 = int(idx3)
        tri_list.append([idx1,idx2,idx3])
        
    #################
    #this generates a fillable collection of tesselations
    x_list = []
    y_list = []   
    
    for t in tri_list:
        if all([ cpt[t[i]]<threshold for i in range(3)]):
            short_x_list = [ xpt[t[i]] for i in range(3)]
            short_y_list = [ ypt[t[i]] for i in range(3)]
            x_list.extend(short_x_list)
            x_list.append(xpt[t[0]])
            x_list.append(None)
            y_list.extend(short_y_list)
            y_list.append(ypt[t[0]])
            y_list.append(None)
    
    return (x_list,y_list)


if __name__=='__main__':
    
    import random
    
    N=100
    xpt = [random.random()-0.5 for i in range(0,N)]
    ypt = [random.random()-0.5 for i in range(0,N)]
    cpt = [random.random()-0.5 for i in range(0,N)]
    
    import matplotlib.pyplot as pp
    
    pp.figure()
    (x_list,y_list) = delaunay2D(xpt,ypt,cpt)
    pp.plot(xpt,ypt,'b.')
    pp.plot(x_list,y_list,'k-')

    pp.fill(x_list,y_list,'g',alpha=0.25,edgecolor='none')    

    pp.show()
    

---

When run, this module will produce a diagram with the cells with filled triangulation where the vertices have a random value greater than 0.

Test Conditions

• PythonXY 2.5.6.2
• Matplotlib 1.0.0
• Qhull 2010.1

---

Answers

• How to tesselate a set of points
• How to generate a Delaunay Triangulation
• How to use Qhull with Python

This work is licensed under a Creative Commons Attribution By license.

Rough Draft: How to Generate Voronoi Diagrams in Python & Matplotlib

Voronoi diagrams partition a space into cells. Each cell represents the region nearest a point  in a set. Voronoi diagrams can be useful for visualizing spatial data.

Qhull is a program which can generate tesselations, convex hulls, and vonoroi diagrams from a set of points. This program is available as a precompiled executable and source code. By interfacing to the command line version of this program, a Voronoi diagram can be generated in Matplotlib.

Sample CodE

To use this code, download Qhull and copy the ‘qvoronoi.exe’ to the working directory. Make sure the code has read, write, and execute privileges in the working directory. Make sure there are not files named ‘data.txt’ or ‘results.txt’ which need to be preserved.

'''
This module generates 2D voronoi diagrams from a list of points.
The vornoi cells can be filtered by supplying a value associated with 
   each node.
'''

def voronoi2D(xpt,ypt,cpt=None,threshold=0):
    '''
    This function returns a list of line segments which describe the voronoi
        cells formed by the points in zip(xpt,ypt).
    
    If cpt is provided, it identifies which cells should be returned.
        The boundary of the cell about (xpt[i],ypt[i]) is returned 
            if cpt[i]<=threshold.
            
    This function requires qvoronoi.exe in the working directory. 
    The working directory must have permissions for read and write access.
    This function will leave 2 files in the working directory:
        data.txt
        results.txt
    This function will overwrite these files if they already exist.
    '''
    
    if cpt is None:
        # assign a value to cpt for later use
        cpt = [0 for x in xpt]
    
    # write the data file
    pts_filename = 'data.txt'
    pts_F = open(pts_filename,'w')
    pts_F.write('2 # this is a 2-D input set\n')
    pts_F.write('%i # number of points\n' % len(xpt))
    for i,(x,y) in enumerate(zip(xpt,ypt)):
        pts_F.write('%f %f # data point %i\n' % (x,y,i))
    pts_F.close()

    # trigger the shell command
    import subprocess
    p = subprocess.Popen('qvoronoi TI data.txt TO results.txt p FN Fv QJ', shell=True)
    p.wait()

    # open the results file and parse results
    results = open('results.txt','r')

    # get 'p' results - the vertices of the voronoi diagram
    data = results.readline()
    voronoi_x_list = []
    voronoi_y_list = []
    data = results.readline()
    for i in range(0,int(data)):
        data = results.readline()
        xx,yy,dummy = data.split(' ')    
        voronoi_x_list.append(float(xx))
        voronoi_y_list.append(float(yy))
        
    # get 'FN' results - the voronoi edges
    data = results.readline()
    voronoi_idx_list = []
    for i in range(0,int(data)):
        data = results.readline()
        this_list = data.split(' ')[:-1]
        for j in range(len(this_list)):
            this_list[j]=int(this_list[j])-1
        voronoi_idx_list.append(this_list[1:])
        
    # get 'FV' results - pairs of points which define a voronoi edge
    # combine these results to build a complete representation of the 
    data = results.readline()
    voronoi_dict = {}
    for i in range(0,int(data)):
        data = results.readline().split(' ')

        pair_idx_1 = int(data[1])
        pair_idx_2 = int(data[2])

        vertex_idx_1 = int(data[3])-1
        vertex_idx_2 = int(data[4])-1

        try:
            voronoi_dict[pair_idx_1].append({ 'edge_vertices':[vertex_idx_1,vertex_idx_2],
                                      'neighbor': pair_idx_2 })
        except KeyError:
            voronoi_dict[pair_idx_1] = [{ 'edge_vertices':[vertex_idx_1,vertex_idx_2],
                                      'neighbor': pair_idx_2 } ]

        try:
            voronoi_dict[pair_idx_2].append({ 'edge_vertices':[vertex_idx_1,vertex_idx_2],
                                      'neighbor': pair_idx_1 })
        except KeyError:
            voronoi_dict[pair_idx_2] = [{ 'edge_vertices':[vertex_idx_1,vertex_idx_2],
                                      'neighbor': pair_idx_1 } ]    

                    
    #################
    # generate a collection of voronoi cells
    x_list = []
    y_list = []    
    for point_idx in voronoi_dict.keys():
        if cpt[point_idx]<=threshold:
            # display this cell, so add the data to the edge list
            e_list = []
            for edge in voronoi_dict[point_idx]:
                p1_idx = edge['edge_vertices'][0]
                p2_idx = edge['edge_vertices'][1]
                e_list.append((p1_idx,p2_idx))
            
            # put the vertices points in order so they
            #   walk around the voronoi cells
            p_list = [p1_idx]
            while True:
                p=p_list[-1]
                for e in e_list:
                    if p==e[0]:
                        next_p = e[1]
                        break
                    elif p==e[1]:
                        next_p = e[0]
                        break
                p_list.append(next_p)
                e_list.remove(e)
                if p_list[0]==p_list[-1]:
                    # the cell is closed
                    break
                
            # build point list
            if all([p>=0 for p in p_list]):
                for p in p_list:
                    if p>=0:
                        x_list.append(voronoi_x_list[p])
                        y_list.append(voronoi_y_list[p])                    
            x_list.append(None)
            y_list.append(None)
                    
    return (x_list,y_list)

if __name__=='__main__':
    
    import random
    
    xpt = [random.random()-0.5 for i in range(0,100)]
    ypt = [random.random()-0.5 for i in range(0,100)]
    cpt = [random.random()-0.5 for i in range(0,100)]
    
    (x_list,y_list) = voronoi2D(xpt,ypt,cpt)
    
    import matplotlib.pyplot as pp
    
    pp.figure()
    pp.plot(xpt,ypt,'b.')
    pp.plot(x_list,y_list,'k-')
    
    pp.fill(x_list,y_list,'g',alpha=0.25,edgecolor='none')    
     
    pp.axis([-1,1,-1,1])

    pp.show()
    

---

When run, this module will produce a Voronoi diagram with the cells which have a value greater than zero visible.

Test Conditions

• PythonXY 2.5.6.2
• Matplotlib 1.0.0
• Qhull 2010.1

This work is licensed under a Creative Commons Attribution By license.

How to quickly plot polygons in Matplotlib

Problem statement:

You need to plot a large collection of polygons in Matplotlib.

Solution:

A simple way to plot filled polygons in Matplotlib is to use the fill function. If you try to plot a collection of  polygons in Matplotlib using sequential calls to fill, it can take a lot of time to generate the graph. There are two ways to speed up the plotting. The first is to use Python’s extended call syntax and pass multiple polygons at one time. The other is to create a single list of points which are the corners of the polygon and separate each polygon by a ‘None’ entry. The extended call method yields an appreciable improvement in performance. Forming a single list of polygon corners separated by ‘None’ results in a significant time savings. For the test cases below, using extended call syntax decreased drawing time by about 50% and using a single list separated by ‘None’ decreased drawing time by about 99%.

Execution time for sequential plotting = 2.496564 sec
Execution time for extended call plotting = 1.172675 sec  
Execution time when using None = = 0.028234 sec

Sample Code:

'''
Code snippet to plot multiple triangle's in Matplotlib using different methods.
'''
import time
import sys
import matplotlib.pyplot as pp
import random

if sys.platform == "win32":
     # On Windows, the best timer is time.clock()
     default_timer = time.clock 
else:     
    # On most other platforms the best timer is time.time()     
    default_timer = time.time
    
    
# generate ends for the triangle line segments
xtris = []
ytris = []
for i in range(1000):
    x1 = random.random()
    x2 = x1 + random.random()/5.0
    x3 = x1 + random.random()/5.0
    xtips = [x1,x2,x3]
    y1 = random.random()
    y2 = y1 + random.random()/5.0
    y3 = y1 + random.random()/5.0
    ytips = [y1,y2,y3]
    xtris.append(xtips)
    ytris.append(ytips)

############################
# time sequential call to matplotlib
pp.figure()
pp.subplot(1,3,1)

t0 = default_timer()
for xtips,ytips in zip(xtris,ytris):
    pp.fill(xtips,ytips,
            facecolor='b',alpha=0.1, edgecolor='none')
t1 = default_timer()

pp.title('Sequential Plotting')

print 'Execution time for sequential plotting = %f sec' % (t1-t0)

# rebuild ends using none to separate polygons
xlist = []
ylist = []
for xtips,ytips in zip(xtris,ytris):
    xlist.extend(xtips)
    xlist.append(None)
    ylist.extend(ytips)
    ylist.append(None)

############################
# build argument list
call_list = []
for xtips,ytips in zip(xtris,ytris):
    call_list.append(xtips)
    call_list.append(ytips)
    call_list.append('-b')
    
############################
# time single call to matplotlib
pp.subplot(1,3,2)

t0 = default_timer()
pp.fill(*call_list,
            facecolor='b',alpha=0.1, edgecolor='none')

t1 = default_timer()

pp.title('Single Plot extended call')

print 'Execution time for extended call plotting = %f sec' % (t1-t0)

    
############################
# time single call to matplotlib
pp.subplot(1,3,3)

t0 = default_timer()
pp.fill(xlist,ylist,
        facecolor='b',alpha=0.1,edgecolor='none')
t1 = default_timer()

pp.title('Single Plot Using None')

print 'Execution time when using None = %f sec' % (t1-t0)

pp.show()

---

Discussion:

   Sequential call and extended call syntax produce the same plot results. However, the single list of points using ‘None’ can generate a different image if transparency is used. In the sample code, ‘alpha=0.1’ was used. When sequential fill or extended call syntax is used, the color from each polygon is additive. When the polygon corners are combined into a single list of points, the polygon colors are not additive. This is only an issue when using transparency. If alpha is set to 1.0 (or not set at all), then all of the plots will be identical.

---

Test Conditions:

• Win 7, 64 bit
• Intel i7
• PythonXY 2.5.6.2
  • Matplotlib 1.0.0

---

Answers:

• How to speed up Matplotlib plots
• How to plot polygons in Matplotlib

This work is licensed under a Creative Commons Attribution By license.