python - 从椭圆生成数组

Question

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

到目前为止，我已经尝试创建一个空数组，其大小我希望通过所有 x,y 位置计算 x^2/a^2 + y^2/b^2 = 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() 那样需要更大的二维中间阵列。

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

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

虽然 ellipse 是一个 bool 数组(在椭圆内/椭圆外)，但在计算中 False 被视为 0，True 被视为 1(例如 ellipse*10 产生一个包含 0 和 10 值的数组)，因此您可以在你的 NumPy 计算。