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.sstatic.net/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.