# Convolute a gaussian kernel with a large array of off-grid centroids without looping? (how to make “A Thousand (Gaussian) Points of Light” )

For a finite object size diffraction simulator, I need to generate arrays which are the sum of thousands of instances of a Gaussian (or other) 2D kernel at centroids that will not fall in any repeatable way with respect to the grid points.

Below is a simple example with a simple hexagonal arrangement of centroids for clarity, along with a much faster analytical expression that comes close to approximating this particularly simple example, but in general the arrangement will be more complicated or even random.

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()


If you simply want to compute the convolution you can construct a Kernel matrix for computation with arbitrary arrangements of Gaussians that do not lie on grid points of the domain where you want to compute the convolution.

To do that one can employ an "outer sum". In python its given as done here: https://stackoverflow.com/questions/33848599/performing-outer-addition-with-numpy

In essense, one can discretize a Fredholm integral equation of the first kind into a system of linear equations to describe the convolution, i.e. in the form of $$Aw=b$$ and the matrix $$A$$ can be computed via feeding this outer sum $$r_c−r_g^T$$ into the Gaussian. The vector $$r_c$$ has all the positions in the domain you want to compute the convolution and $$r_g$$ is the vector that containts all central positions of the gaussians. The weight vector $$w$$ describes the amplitudes of the gaussians and $$b$$ is the resultant convolution you seek.

I rewrote your example this way in Matlab below:

%%% 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)

• This is really great, thank you! I'll give the method a try and also do some reading. I think that this is going to be one of those answers that keeps on giving... – uhoh Jul 21 at 13:26
• – uhoh Jul 21 at 13:28
• You're welcome. This approach also solves your other question with the parametric function. In this case, the weight vector are the function values in $g(t)$ – Ron Jul 21 at 13:29