我目前正在尝试编写自己的 2D 高斯函数作为编码练习，并且已经能够创建以下脚本：
import numpy as np
import matplotlib.pyplot as plt
def Gaussian2D_v1(coords=None, # x and y coordinates for each image.
amplitude=1, # Highest intensity in image.
xo=0, # x-coordinate of peak centre.
yo=0, # y-coordinate of peak centre.
sigma_x=1, # Standard deviation in x.
sigma_y=1, # Standard deviation in y.
rho=0, # Correlation coefficient.
offset=0): # Offset from zero (background radiation).
x, y = coords
xo = float(xo)
yo = float(yo)
# Create covariance matrix
mat_cov = [[sigma_x**2, rho * sigma_x * sigma_y],
[rho * sigma_x * sigma_y, sigma_y**2]]
mat_cov = np.asarray(mat_cov)
# Find its inverse
mat_cov_inv = np.linalg.inv(mat_cov)
G_array = []
# Calculate pixel by pixel
# Iterate through row last
for i in range(0, np.shape(y)[0]):
# Iterate through column first
for j in range(0, np.shape(x)[1]):
mat_coords = np.asarray([[x[i, j]-xo],
[y[i, j]-xo]])
G = (amplitude * np.exp(-0.5*np.matmul(np.matmul(mat_coords.T,
mat_cov_inv),
mat_coords)) + offset)
G_array.append(G)
G_array = np.asarray(G_array)
G_array = G_array.reshape(64, 64)
return G_array.ravel()
coords = np.meshgrid(np.arange(0, 64), np.arange(0, 64))
model_1 = Gaussian2D_v1(coords,
amplitude=20,
xo=32,
yo=32,
sigma_x=6,
sigma_y=3,
rho=0.8,
offset=20).reshape(64, 64)
plt.figure(figsize=(5, 5)).add_axes([0,
0,
1,
1])
plt.contourf(model_1)
代码是有效的，但如您所见，我目前正在一次一个点地遍历网格，并将每个点附加到列表中，然后将其转换为数组并重新调整形状以给出 2D 高斯分布。
如何修改脚本以放弃使用嵌套的"for"循环，并让程序考虑整个网格进行矩阵计算？这种方法可能吗？
谢谢！
当然有一个解决方案，numpy 是关于数组操作和代码的矢量化！np.matmul可以采用超过 2 个维度的参数，并仅在最后两个轴上应用矩阵乘法(并且此计算在其他轴上并行(。但是，确保正确的轴顺序可能会变得棘手。
这是您编辑的代码：
import numpy as np
import matplotlib.pyplot as plt
def Gaussian2D_v1(coords, # x and y coordinates for each image.
amplitude=1, # Highest intensity in image.
xo=0, # x-coordinate of peak centre.
yo=0, # y-coordinate of peak centre.
sigma_x=1, # Standard deviation in x.
sigma_y=1, # Standard deviation in y.
rho=0, # Correlation coefficient.
offset=0): # Offset from zero (background radiation).
x, y = coords
xo = float(xo)
yo = float(yo)
# Create covariance matrix
mat_cov = [[sigma_x**2, rho * sigma_x * sigma_y],
[rho * sigma_x * sigma_y, sigma_y**2]]
mat_cov = np.asarray(mat_cov)
# Find its inverse
mat_cov_inv = np.linalg.inv(mat_cov)
# PB We stack the coordinates along the last axis
mat_coords = np.stack((x - xo, y - yo), axis=-1)
G = amplitude * np.exp(-0.5*np.matmul(np.matmul(mat_coords[:, :, np.newaxis, :],
mat_cov_inv),
mat_coords[..., np.newaxis])) + offset
return G.squeeze()
coords = np.meshgrid(np.arange(0, 64), np.arange(0, 64))
model_1 = Gaussian2D_v1(coords,
amplitude=20,
xo=32,
yo=32,
sigma_x=6,
sigma_y=3,
rho=0.8,
offset=20)
plt.figure(figsize=(5, 5)).add_axes([0, 0, 1, 1])
plt.contourf(model_1)
因此，方程是exp(-0.5 * (X - μ(' Cinv (X - μ( (，其中X是我们的坐标矩阵，μ平均值(x0，y0(，Cinv是逆协方差矩阵('是转置(。在代码中，我将两个网格堆叠到一个新的矩阵中，以便：mat_coords的形状为 (Ny， Nx， 2(。在第一个np.matmul调用中，我添加了一个新轴，以便形状类似于:(Ny， Nx， 1， 2( * (2， 2( = (Ny， Nx， 1， 2(。如您所见，矩阵乘法是在最后两个轴上完成的，在另一个轴上并行完成。然后，我添加一个新轴，以便：(Ny， Nx， 1， 2( * (Ny， Nx， 2， 1( = (Ny， Nx， 1， 1(。np.squeeze()调用返回一个没有最后两个单例轴的版本。