从椭圆生成数组

Question

我有一个方程，它以一般形式 x^2/a^2 + y^2/b^2 = 1 创建一个椭圆。我希望生成一个数组，其中椭圆内的所有点都设置为 1，并且外面的所有点都是零。然后将这个数组与另一个数组进行卷积。

到目前为止，我已经尝试创建一个空数组，其大小我希望遍历所有 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) 为中心，因此我希望它们位于数组的中心，但它们出现在左上角。

有没有人知道我哪里出错了？或者也许是生产我的阵列的更好方法？

---编辑---

尝试过 EOL 的回答后，我收到了一个椭圆数组，但它与它应该建模的椭圆不匹配。这是一张图片来说明我的意思：http://i.imgur.com/M4vh4il.jpg?1 椭圆数组没有椭圆的旋转。 生成椭圆和椭圆数组的代码如下：

    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，而不是您在OP 中使用的像素数。因为我的坐标系以 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() 更好。

1234563 987654326@ 仍然是行向量）。也不需要像 numpy.meshgrid() 那样需要更大的 2D 中间数组。

• NumPy 可以执行“在椭圆中吗？”直接测试（即 <= 1 部分），没有 Python 循环。

现在，如果你只能使用整数坐标，x 和 y 可以简单地使用 np.arange() 而不是更“精致”的np.linspace() 获得。

虽然ellipse 是一个布尔数组（椭圆内/外），但在计算中，False 被视为 0，True 被视为 1（例如，ellipse*10 产生一个包含 0 和 10 值的数组），所以你可以在你的 NumPy 计算中使用它。