Extract coordinates enclosed by a matplotlib patch.

Question

I have created an ellipse using matplotlib.patches.ellipse as shown below:

patch = mpatches.Ellipse(center, major_ax, minor_ax, angle_deg, fc='none', ls='solid', ec='g', lw='3.')

What I want is a list of all the integer coordinates enclosed inside this patch. I.e. If I was to plot this ellipse along with every integer point on the same grid, how many of those points are enclosed in the ellipse?

I have tried seeing if I can extract the equation of the ellipse so I can loop through each point and see whether it falls within the line but I can't seem to find an obvious way to do this, it becomes more complicated as the major axis of the ellipse can be orientated at any angle. The information to do this must be stored in patches somewhere, but I can't seem to find it.

Any advice on this would be much appreciated.

Answer 1

Ellipse objects have a method contains_point which will return 1 if the point is in the ellipse, 0 other wise.

Stealing from @DrV 's answer:

import matplotlib.pyplot as plt
import matplotlib.patches
import numpy as np

# create an ellipse
el = matplotlib.patches.Ellipse((50,-23), 10, 13.7, 30, facecolor=(1,0,0,.2), edgecolor='none')

# calculate the x and y points possibly within the ellipse
y_int = np.arange(-30, -15)
x_int = np.arange(40, 60)

# create a list of possible coordinates
g = np.meshgrid(x_int, y_int)
coords = list(zip(*(c.flat for c in g)))

# create the list of valid coordinates (from untransformed)
ellipsepoints = np.vstack([p for p in coords if el.contains_point(p, radius=0)])

# just to see if this works
fig = plt.figure()
ax = fig.add_subplot(111)
ax.add_artist(el)
ep = np.array(ellipsepoints)
ax.plot(ellipsepoints[:,0], ellipsepoints[:,1], 'ko')
plt.show()

This will give you the result as below:

Answer 2

If you really want to use the methods offered by matplotlib, then:

import matplotlib.pyplot as plt
import matplotlib.patches
import numpy as np

# create an ellipse
el = matplotlib.patches.Ellipse((50,-23), 10, 13.7, 30, facecolor=(1,0,0,.2), edgecolor='none')

# find the bounding box of the ellipse
bb = el.get_window_extent()

# calculate the x and y points possibly within the ellipse
x_int = np.arange(np.ceil(bb.x0), np.floor(bb.x1) + 1, dtype='int')
y_int = np.arange(np.ceil(bb.y0), np.floor(bb.y1) + 1, dtype='int')

# create a list of possible coordinates
g = np.meshgrid(x_int, y_int)
coords = np.array(zip(*(c.flat for c in g)))

# create a list of transformed points (transformed so that the ellipse is a unit circle)
transcoords = el.get_transform().inverted().transform(coords)

# find the transformed coordinates which are within a unit circle
validcoords = transcoords[:,0]**2 + transcoords[:,1]**2 < 1.0

# create the list of valid coordinates (from untransformed)
ellipsepoints = coords[validcoords]

# just to see if this works
fig = plt.figure()
ax = fig.add_subplot(111)
ax.add_artist(el)
ep = np.array(ellipsepoints)
ax.plot(ellipsepoints[:,0], ellipsepoints[:,1], 'ko')

Seems to work:

(Zooming in reveals that even the points hanging on the edge are inside.)

The point here is that matplotlib handles ellipses as transformed circles (translate, rotate, scale, anything affine). If the transform is applied in reverse, the result is a unit circle at origin, and it is very simple to check if a point is within that.

Just a word of warning: The get_window_extent may not be extremely reliable, as it seems to use the spline approximation of a circle. Also, see tcaswell's comment on the renderer-dependency.

In order to find a more reliable bounding box, you may:

• create a horizontal and vertical vector into the plot coordinates (their position is not important, ([0,0],[1,0]) and ([0,0], [0,1]) will do)

• transform these vectors into the ellipse coordinates (the get_transform, etc.)

• find in the ellipse coordinate system (i.e. the system where the ellipse is a unit circle around the origin) the four tangents of the circle which are parallel to these two vectors

• find the intersection points of the vectors (4 intersections, but 2 diagonal will be enough)

• transform the intersection points back to the plot coordinates

This will give an accurate (but of course limited by the numerical precision) square bounding box.

However, you may use a simple approximation:

• all possible points are within a circle whose center is the same as that of the ellipse and whose diameter is the same as that of the major axis of the ellipse

In other words, all possible points are within a square bounding box which is between x0+-m/2, y0+-m/2, where (x0, y0) is the center of the ellipse and m the major axis.

Answer 3

I'd like to offer another solution that uses the Path object's contains_points() method instead of contains_point():

First get the coordinates of the ellipse and make it into a Path object:

elpath=Path(el.get_verts())

(NOTE that el.get_paths() won't work for some reason.)

Then call the path's contains_points():

validcoords=elpath.contains_points(coords)

Below I'm comparing @tacaswell's solution (method 1), @Drv's (method 2) and my own (method 3) (I've enlarged the ellipse by ~5 times):

import numpy
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from matplotlib.path import Path
import time

#----------------Create an ellipse----------------
el=Ellipse((50,-23),50,70,30,facecolor=(1,0,0,.2), edgecolor='none')

#---------------------Method 1---------------------
t1=time.time()

for ii in range(50):
    y=numpy.arange(-100,50)
    x=numpy.arange(-30,130)
    g=numpy.meshgrid(x,y)
    coords=numpy.array(zip(*(c.flat for c in g)))

    ellipsepoints = numpy.vstack([p for p in coords if el.contains_point(p, radius=0)])

t2=time.time()
print 'time of method 1',t2-t1

#---------------------Method 2---------------------
t2=time.time()

for ii in range(50):
    y=numpy.arange(-100,50)
    x=numpy.arange(-30,130)
    g=numpy.meshgrid(x,y)
    coords=numpy.array(zip(*(c.flat for c in g)))

    invtrans=el.get_transform().inverted()
    transcoords=invtrans.transform(coords)
    validcoords=transcoords[:,0]**2+transcoords[:,1]**2<=1.0
    ellipsepoints=coords[validcoords]

t3=time.time()
print 'time of method 2',t3-t2

#---------------------Method 3---------------------
t3=time.time()

for ii in range(50):
    y=numpy.arange(-100,50)
    x=numpy.arange(-30,130)
    g=numpy.meshgrid(x,y)
    coords=numpy.array(zip(*(c.flat for c in g)))

    #------Create a path from ellipse's vertices------
    elpath=Path(el.get_verts())
    # call contains_points()
    validcoords=elpath.contains_points(coords)
    ellipsepoints=coords[validcoords]

t4=time.time()
print 'time of method 3',t4-t3

#---------------------Plot it ---------------------
fig,ax=plt.subplots()
ax.add_artist(el)
ep=numpy.array(ellipsepoints)
ax.plot(ellipsepoints[:,0],ellipsepoints[:,1],'ko')
plt.show(block=False)

I got these execution time:

time of method 1 62.2502269745
time of method 2 0.488734006882
time of method 3 0.588987112045

So the contains_point() approach is way slower. The coordinate-transformation method is faster than mine, but when you get irregular shaped contours/polygons, this method would still work.

Finally the result plot: