计算科学 - 用大量离网质心卷积一个高斯核而不循环？（如何制作“一千（高斯）光点”） - 吾爱随笔录

用大量离网质心卷积一个高斯核而不循环？（如何制作“一千（高斯）光点”）

计算科学 Python 网格生成卷积

2021-12-19 22:06:09

对于有限物体大小的衍射模拟器，我需要生成数组，这些数组是质心处高斯（或其他）2D 内核的数千个实例的总和，这些实例不会以任何可重复的方式落在网格点上。

下面是一个简单的例子，为了清楚起见，质心简单的六边形排列，以及一个更快的解析表达式，接近于这个特别简单的例子，但通常排列会更复杂甚至是随机的。

import numpy as np
import matplotlib.pyplot as plt

hw = 10
N = 250
twoNp1 = 2 * N + 1
ximg, yimg = np.mgrid[-hw:hw:twoNp1 * 1j, -hw:hw:twoNp1 * 1j]
a = np.sqrt(2)
sig1 = hw/3.
sig2 = a/5.

r3o2, twopi = np.sqrt(3) / 2, 2 * np.pi

vecs = a * np.array([[1, 0], [1/2, r3o2]])
nmax = int(2 * (ximg.max()/a)) # overfill overkill
i, j = [thing.flatten() for thing in np.mgrid[-nmax:nmax+1, -nmax:nmax+1]]
keep = np.abs(i + j) <= nmax
i, j = [thing[keep, None] for thing in (i, j)]  
points = i * vecs[0] + j * vecs[1]

img = np.zeros_like(ximg)
for x, y in points:
    img += np.exp(-((x-ximg)**2 + (y-yimg)**2) / (2 * sig2**2))

# actual gaussians
img *= np.exp(-((ximg)**2 + (yimg)**2) / (2 * sig1**2))

# sinusoidal pattern
k = (twopi / (a * r3o2)) * np.array([[1, 0], [0.5, r3o2], [0.5, -r3o2]])
# that's 0, 60, -60
# try -30, 30, 90
k = (twopi / (a * r3o2)) * np.array([[r3o2, -1/2], [r3o2, 1/2], [0, 1]])

phases = [kay[0] * ximg + kay[1] * yimg for kay in k]
amplitudes = [np.cos(phase) for phase in phases]
amplitude = (1.5 + sum(amplitudes))/4.5
amplitude *= np.exp(-((ximg)**2 + (yimg)**2) / (2 * sig1**2))

if False:
    n = int((2 * hw / a)**2 + 0.5)
    randoms = hw * (2 * np.random.random((n,2)) - 1)
    rand = np.zeros_like(ximg)
    for x, y in randoms:
        rand += np.exp(-((x-ximg)**2 + (y-yimg)**2) / (2 * sig2**2))
    rand *= np.exp(-((ximg)**2 + (yimg)**2) / (2 * sig1**2))

if True:
    fig, (ax1, ax2) = plt.subplots(2, 1)
    extent = [ximg.min(), ximg.max(), yimg.min(), yimg.max()]
    one = ax1.imshow(img, origin='lower', extent=extent)
    ax1.plot(ximg, hw * (img[N] - 1), '-r', linewidth=0.5)
    ax1.set_title('"thousands" of Gaussians')
    # fig.colorbar(one, ax=ax1)
    two = ax2.imshow(amplitude, origin='lower', extent=extent)
    ax2.plot(ximg, hw * (amplitude[N] - 1), '-r', linewidth=0.5)
    ax2.set_title('sinusoidal approximation')
    # fig.colorbar(two, ax=ax2)
    plt.show()

3个回答

如果您只是想计算卷积，您可以构建一个核矩阵，用于计算任意排列的高斯矩阵，这些高斯矩阵不位于要计算卷积的域的网格点上。

为此，可以使用“外部总和”。在 python 中，它在这里完成： https ://stackoverflow.com/questions/33848599/performing-outer-addition-with-numpy

本质上，可以将第一类 Fredholm 积分方程离散化为线性方程组来描述卷积，即形式为 $Aw=b$ 和矩阵 $A$ 可以通过输入这个外部总和来计算 $r_c−r_g^T$ 进入高斯。向量 $r_c$ 具有要计算卷积的域中的所有位置，并且 $r_g$ 是包含所有高斯中心位置的向量。权重向量 $w$ 描述了高斯的幅度和 $b$ 是您寻求的结果卷积。

我在下面的 Matlab 中以这种方式重写了您的示例：

%%% Positions to compute convolution
N = 150;
x = linspace(-10,10,N); 
y = linspace(-10,10,N);
[X,Y] = ndgrid(x,y);

%%% Gaussian Positions
B = [1,cosd(120);0,sind(120)];
m = 30;
d = 1.5;
[i,j] = ndgrid(-m:d:m,-m:d:m);
P = B*[i(:)';j(:)'];

%%% Compute Kernel A
Xb = X(:) - P(1,:);
Yb = Y(:) - P(2,:);
sig1 = 5;
A = exp(-(Xb.^2+Yb.^2)*sig1); 

%%% Compute Convolution
sig2 = 2e-2;
w = ones(size(A,2),1); % weight of individual gaussians
b = reshape(A*w,size(X));
% b = b.*exp(-(X.^2+Y.^2)*sig2) % if needed

%%% Plotting
imagesc(b)

大约在 Python 中实现@Ron 的答案，关键是img = (A * weights).sum(axis=2)哪个当然是循环，但它是在 numpy 调用的编译代码中完成的。

这里可以通过调整数组定义来优化速度，以便我们可以对不同的轴和其他东西进行求和。查看 Code Review 的答案

出于演示目的，我添加了随机权重。

用大量离网质心卷积一个高斯核而不循环？

可能是一个有点误导性的标题，因为问题实际上是有限数量的单个质心的总和。相关的（目前尚未回答的）问题将更具挑战性，因为它是一个真正的卷积。

如何在 Python 中生成 f(x, y) 与参数函数 g(t), x(t), y(t) 的卷积？（比这个蛮力总和更好的东西）

import numpy as np
import matplotlib.pyplot as plt

# Positions to compute convolution
N = 150
x = np.linspace(-10, 10, 150)
y = np.linspace(-10, 10, 150)
X, Y = np.meshgrid(x, y)

# Gaussian Positions
v1, v2 = np.array([[1, 0], [0.5,np.sqrt(3)/2]])
m = 30
i = np.linspace(-m, m, 40) 
I, J = [thing.flatten() for thing in np.meshgrid(i, i)]
P = v1[:, None] * I + v2[:, None] * J
keep = np.abs(I + J) <= m # makes the boundary hexagonal
P = P[:, keep]

Xb = X[..., None] - P[1]
Yb = Y[..., None] - P[0]
sig1 = 5
A = np.exp(-(Xb**2 + Yb**2) * sig1)
weights = np.random.random(len(P[1]))
img = (A * weights).sum(axis=2)

对“源”和“目标”点的任意集合快速执行此操作的经典方法是使用由 Greengard 开发的称为快速高斯变换的快速多极类型算法。一个快速的谷歌搜索出现了这个 Julia 包，实现了作者认为是“改进的”FGT：https ://github.com/jwmerrill/FastGaussTransforms.jl

基于“改进的快速高斯变换”，计算机视觉国际会议论文集，464 -471，C. Yang、R. Duraiswami、NA Gumerov、L. Davis。pdf

其它你可能感兴趣的问题

上一篇通过 SVD 进行线性回归不会产生与递增多项式次数的最佳拟合下一篇实施一个O ( n对数( n ) )O(nlog⁡(n))计算实对称三对角矩阵特征值的方法