Numpy 操作真值数组，取决于 x/y 索引

Question

所以我有一个数组（它很大 - 2048x2048），我想根据它们的位置进行一些元素操作。我很困惑如何做到这一点（我被告知不要使用 for 循环，当我尝试我的 IDE 冻结并且运行速度非常慢时）。

关于问题：

h = aperatureimage
h[:,:] = 0
indices = np.where(aperatureimage>1)

for True in h:
    h[index] = np.exp(1j*k*z)*np.exp(1j*k*(x**2+y**2)/(2*z))/(1j*wave*z)

所以我有一个索引，它（我在这里假设）本质上是我更大的 aperatureimage 数组的“裁剪”版本。 *注意：孔径图像是转换为数组的灰度图像，上面有形状或文字，我想找到孔径的所有“白色”区域并执行我的操作。

如何访问索引的各个 x/y 值，这将允许我执行指数运算？当我尝试 index[:,None] 时，导致程序吐出“ValueError：广播尺寸太大”。我也得到数组不可广播以纠正形状。任何帮助将不胜感激！

再澄清一点：x 和 y 是我想要更改的唯一值（基本上是我的数组中存在白色、z、k 以及之前定义的任何其他值的点）。

编辑：

我不确定我上面发布的代码是否正确，它返回两个空数组。当我这样做时 指数=（孔径图像==1） 打印长度（索引）

实际上，到目前为止，我所做的一切都无法正常工作。我有一个 2048x2048 的图像，中间有一个 128x128 的白色方块。我想将此图像转换为数组，查看所有值并确定数组不是黑色的索引值（x，y）（我只有白色/黑色，双层图像对我不起作用）。然后我想取数组不为 0 的所有值 (x,y)，并将它们乘以上面列出的 h[index] 值。

如有需要，我可以发布更多信息。如果你看不出来，我就卡住了。

EDIT2：这里有一些可能有帮助的代码——我想我已经解决了上面的问题（我现在可以访问数组的成员并对它们执行操作）。但是 - 由于某种原因，我的 for 循环中的 Fx 值永远不会增加，它会永远循环 Fy....

import sys, os
from scipy.signal import *
import numpy as np
import Image, ImageDraw, ImageFont, ImageOps, ImageEnhance, ImageColor



def createImage(aperature, type):
    imsize = aperature*8
    middle = imsize/2
    im = Image.new("L", (imsize,imsize))
    draw = ImageDraw.Draw(im)
    box = ((middle-aperature/2, middle-aperature/2), (middle+aperature/2, middle+aperature/2))
import sys, os
from scipy.signal import *
import numpy as np
import Image, ImageDraw, ImageFont, ImageOps, ImageEnhance, ImageColor



def createImage(aperature, type):

    imsize = aperature*8 #Add 0 padding to make it nice
    middle = imsize/2 # The middle (physical 0) of our image will be the imagesize/2
    im = Image.new("L", (imsize,imsize)) #Make a grayscale image with imsize*imsize pixels
    draw = ImageDraw.Draw(im) #Create a new draw method
    box = ((middle-aperature/2, middle-aperature/2), (middle+aperature/2, middle+aperature/2)) #Bounding box for aperature
    if type == 'Rectangle':
        draw.rectangle(box, fill = 'white') #Draw rectangle in the box and color it white
    del draw
    return im, middle


def Diffraction(aperaturediameter = 1, type = 'Rectangle', z = 2000000, wave = .001):

    # Constants
    deltaF = 1/8 # Image will be 8mm wide
    z = 1/3.
    wave = 0.001
    k = 2*pi/wave

    # Now let's get to work
    aperature = aperaturediameter * 128 # Aperaturediameter (in mm) to some pixels
    im, middle = createImage(aperature, type) #Create an image depending on type of aperature 
    aperaturearray = np.array(im) # Turn image into numpy array

    # Fourier Transform of Aperature
    Ta = np.fft.fftshift(np.fft.fft2(aperaturearray))/(len(aperaturearray)) 

    # Transforming and calculating of Transfer Function Method
    H = aperaturearray.copy() # Copy image so H (transfer function) has the same dimensions as aperaturearray
    H[:,:] = 0 # Set H to 0
    U = aperaturearray.copy()
    U[:,:] = 0
    index = np.nonzero(aperaturearray) # Find nonzero elements of aperaturearray
    H[index[0],index[1]] = np.exp(1j*k*z)*np.exp(-1j*k*wave*z*((index[0]-middle)**2+(index[1]-middle)**2)) # Free space transfer for ap array
    Utfm = abs(np.fft.fftshift(np.fft.ifft2(Ta*H))) # Compute intensity at distance z

    # Fourier Integral Method
    apindex = np.nonzero(aperaturearray)
    U[index[0],index[1]] = aperaturearray[index[0],index[1]] * np.exp(1j*k*((index[0]-middle)**2+(index[1]-middle)**2)/(2*z))
    Ufim = abs(np.fft.fftshift(np.fft.fft2(U))/len(U))

    # Save image
    fim = Image.fromarray(np.uint8(Ufim))
    fim.save("PATH\Fim.jpg")
    ftfm = Image.fromarray(np.uint8(Utfm))
    ftfm.save("PATH\FTFM.jpg")
    print "that may have worked..."
    return

if __name__ == '__main__':
    Diffraction()

您需要 numpy、scipy 和 PIL 才能使用此代码。

当我运行它时，它会遍历代码，但其中没有数据（一切都是黑色的）。现在我在这里遇到了一个真正的问题，因为我不完全理解我正在做的数学（这是针对硬件的），而且我对 Python 没有牢牢掌握。

U[index[0],index[1]] = aperaturearray[index[0],index[1]] * np.exp(1j*k*((index[0]-middle)**2+(index[1]-middle)**2)/(2*z))

那条线应该适用于对我的数组执行元素计算吗？

Answer 1

您能否发布一个最小但完整的示例？一个我们可以复制/粘贴并自己运行的？

与此同时，在您当前示例的前两行中：

h = aperatureimage
h[:,:] = 0

您将“aperatureimage”和“h”都设置为 0。这可能不是您想要的。您可能需要考虑：

h = aperatureimage.copy()

这会生成 aperatureimage 的副本，而您的代码只是将 h 指向与 aperatureimage 相同的数组。所以改变一个会改变另一个。 请注意，复制非常大的数组可能会花费更多内存。

我认为您正在尝试做的是：

import numpy as np

N = 2048
M = 64

a = np.zeros((N, N))
a[N/2-M:N/2+M,N/2-M:N/2+M]=1

x,y = np.meshgrid(np.linspace(0, 1, N), np.linspace(0, 1, N))

b = a.copy()

indices = np.where(a>0)

b[indices] = np.exp(x[indices]**2+y[indices]**2)

或类似的东西。在任何情况下，这都会根据“a”大于 0 的 x/y 坐标在“b”中设置一些值。尝试使用 imshow 将其可视化。祝你好运！

关于编辑

您应该标准化您的输出，使其适合 8 位整数。目前，您的一个数组的最大值远大于 255，而一个数组的最大值要小得多。试试这个：

fim = Image.fromarray(np.uint8(255*Ufim/np.amax(Ufim)))
fim.save("PATH\Fim.jpg")
ftfm = Image.fromarray(np.uint8(255*Utfm/np.amax(Utfm)))
ftfm.save("PATH\FTFM.jpg")

还要考虑 np.zeros_like() 而不是复制和清除 H 和 U。

最后，我个人非常喜欢在开发这样的东西时使用 ipython。如果将代码放在脚本顶层的 Diffraction 函数中（代替 'if __ name __ &c.'），则可以直接从 ipython 访问变量。像 np.amax(Utfm) 这样的快速命令会告诉你确实有值！=0。 imshow() 总是很适合查看矩阵。