Producing an array from an ellipse

Question

I have an equation that creates an ellipse in the general form x^2/a^2 + y^2/b^2 = 1. I wish to produce an array whereby all points inside the ellipse are set to one and all points outside are a zero. This array is then to be convolved with another.

So far I have tried to create an empty array of the size I want the go through all x,y positions calculating x^2/a^2 + y^2/b^2 = 1. If the general form is less than one enter a one into the array, otherwise continue to the next x,y position.

Here is my code:

    arr = numpy.array(im) 
    sh = numpy.shape(arr)
    ar = numpy.empty(sh)

    for x in range (sh[0]):
        xx = x*x
        for y in range (sh[1]):
            yy = y*y
            ellips = xx/(a*a)+yy/(b*b)
            if ellips < 1:
                ar[xx,yy] = '1'
            else:
                break

However, this doesn't produce what I expect as my ellipse is always centered at (0,0) thus I expect the ones to be in the center of my array but they appear in the top left corner.

Does anyone have insight as to where I have gone wrong? Or perhaps a better way to produce my array?

---Edit---

Having tried EOL's answer I receive an ellipse array however it doesn't match the ellipse it should model. Here is a picture to illustrate what I mean: https://i.stack.imgur.com/aVx4I.jpg The ellipse array doesn't have the rotation of the ellipse. The code to produce the ellipse and ellipse array is as follows:

    def Ellipssee(z,stigx,stigy):
        points=100 #Number of points to construct the ellipse
        x0,y0 = 0,0 #Beam is always centred
        z0 = 4 # z0 a constant of the device
        al = 2 # alpha a constant of the device
        de = sqrt((stigx**2 + stigy**2))
        ang = arctan2(stigy, stigx) # result in radians
        a = (z+de-z0)*al
        b = (z-de-z0)*al 
        cos_a,sin_a=cos(ang),sin(ang)
        the=linspace(0,2*pi,points)
        X=a*cos(the)*cos_a-sin_a*b*sin(the)+x0
        Y=a*cos(the)*sin_a+cos_a*b*sin(the)+y0

        img = Image.open("bug.png").convert("L") # load image for array size
        arr = np.array(img) 
        sh = np.shape(arr)

        nx = sh[0]   # number of pixels in x-dir
        ny = sh[1]   # number of pixels in y-dir

        x0 = 0;  # x center, half width                                       
        y0 = 0;  # y center, half height                                      
        x = np.linspace(-60, 60, nx)  # x values of interest
        y = np.linspace(-30, 30, ny)  # y values of interest
        ellipseArr = ((x-x0)/a)**2 + ((y[:,None]-y0)/b)**2 <= 1

I have been calling the method with the values Ellipse(1,6,8).

Why is the rotation being lost in the array creation?

Answer 1

NumPy can do this directly, with no Python loop:

>>> import numpy as np
>>> from matplotlib import pyplot

>>> x0 = 4; a = 5  # x center, half width                                       
>>> y0 = 2; b = 3  # y center, half height                                      
>>> x = np.linspace(-10, 10, 100)  # x values of interest
>>> y = np.linspace(-5, 5, 100)[:,None]  # y values of interest, as a "column" array
>>> ellipse = ((x-x0)/a)**2 + ((y-y0)/b)**2 <= 1  # True for points inside the ellipse

>>> pyplot.imshow(ellipse, extent=(x[0], x[-1], y[0], y[-1]), origin="lower")  # Plot

A few key points of this technique are:

• The squares (in x and y) are calculated only once (and not for each point individually). This is better than using numpy.meshgrid().

• Thanks to NumPy's broadcasting rules, the contributions of x and y are summed together in a simple way (y[:,None] essentially makes y a column vector of y value, while x remain a row vector). There is also no need for larger 2D intermediate array, as would be needed with numpy.meshgrid().

• NumPy can do the "is in ellipse?" test directly (that's the <= 1 part), with no Python loop.

Now, if you can live with integer coordinates only, x and y can simply be obtained with np.arange() instead of the more "refined" np.linspace().

While ellipse is an array of booleans (in/out of the ellipse), False is treated like 0 and True like 1 in calculations (for instance ellipse*10 produces an array of 0 and 10 values), so you can use it in your NumPy calculations.

Answer 2

Here is an answer using some tricks from numpy.

import numpy as np

a = 3.0
b = 1.0
nx = 256   # number of pixels in x-dir
ny = 256   # number of pixels in y-dir

# set up a coordinate system
x = np.linspace(-5.0, 5.0, nx)
y = np.linspace(-5.0, 5.0, ny)

# Setup arrays which just list the x and y coordinates
xgrid, ygrid = np.meshgrid(x, y)

# Calculate the ellipse values all at once
ellipse = xgrid**2 / a**2 + ygrid**2 / b**2

# Create an array of int32 zeros
grey = np.zeros((nx,ny), dtype=np.int32)

# Put 1's where ellipse is less than 1.0
# Note ellipse <1.0 produces a boolean array that then indexes grey
grey[ellipse < 1.0] = 1

Note that I've used coordinates xgrid and ygrid as opposed to pixel numbers as you did in your OP. Because my coordinate system is centered on 0.0,0.0, I get an ellipse in the middle of my array. In your question you used pixel indices which have 0,0 in the corner so you get an ellipse in the corner. You can either move your coordinate system (as I did) or account for the shift in your ellipse equation (as the comment suggested).