从椭圆生成数组

Question

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

到目前为止，我已经尝试创建一个大小的空数组，我希望通过所有x，y位置计算x ^ 2 / a ^ 2 + y ^ 2 / b ^ 2 = 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可以直接执行此操作，没有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成为y值的列向量，而x仍然是行向量）。如numpy.meshgrid()所需，也不需要更大的2D中间阵列。
• NumPy可以做“椭圆形吗？”直接测试（这是<= 1部分），没有Python循环。

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

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

Answer 2

以下是使用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的像素索引，所以你在角落里得到一个椭圆。您可以移动坐标系（就像我一样）或考虑椭圆方程的移位（如评论所示）。