从椭圆生成数组

Question

我有一个公式可以画出一个椭圆，公式是 x^2/a^2 + y^2/b^2 = 1。我想生成一个数组，数组里所有在椭圆内部的点都设置为1，而在外部的点则设置为0。然后这个数组还要和另一个数组进行卷积运算。

到目前为止，我尝试创建一个空数组，大小是我想要的，然后遍历所有的 x 和 y 位置，计算 x^2/a^2 + y^2/b^2 = 1。如果这个公式的结果小于1，就在数组里填入1，否则就继续下一个 x 和 y 的位置。

这是我的代码：

    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

但是，这样做并没有得到我期待的结果，因为我的椭圆总是以 (0,0) 为中心，所以我希望数组的中心位置是1，但实际上它们出现在了左上角。

有没有人能告诉我我哪里出错了？或者有没有更好的方法来生成我的数组？

---编辑---

我尝试了 EOL 的回答，得到了一个椭圆数组，但它并没有和我想要的椭圆匹配。这里有一张图片来说明我的意思： https://i.stack.imgur.com/aVx4I.jpg。这个椭圆数组没有椭圆的旋转效果。

生成椭圆和椭圆数组的代码如下：

    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

我一直用 Ellipse(1,6,8) 这个方法来调用。

为什么在创建数组的时候旋转效果丢失了呢？

Answer 1

这里有一个使用了 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

注意，我使用了坐标 xgrid 和 ygrid，而不是你在问题中提到的像素编号。因为我的坐标系统是以 0.0,0.0 为中心的，所以我在数组的中间得到了一个椭圆。而你在提问时使用的是像素索引，0,0 在角落，所以你得到的椭圆就在角落里。你可以选择移动你的坐标系统（就像我做的那样），或者在你的椭圆方程中考虑这个偏移（正如评论所建议的那样）。

Answer 2

NumPy可以直接做到这一点，不需要用Python的循环：

>>> 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

这个方法有几个关键点：

• x和y的平方只计算一次（而不是每个点单独计算）。这比使用numpy.meshgrid()要好。

• 得益于NumPy的广播规则，x和y的值可以简单地加在一起（y[:,None]实际上把y变成了一个列向量，而x保持为行向量）。这样就不需要像numpy.meshgrid()那样生成更大的二维中间数组。

• NumPy可以直接进行“是否在椭圆内？”的测试（就是<= 1的部分），同样不需要Python的循环。

现在，如果你只需要整数坐标，x和y可以直接用np.arange()来获取，而不是使用更“精细”的np.linspace()。

而ellipse是一个布尔数组（表示在椭圆内或外），在计算中，False被当作0，True被当作1（例如ellipse*10会产生一个包含0和10的数组），所以你可以在NumPy的计算中使用它。