# 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 + j * vecs

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 * ximg + kay * 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
– uhoh
Aug 15 at 4:39

Approximately implementing @Ron's answer in Python the key is img = (A * weights).sum(axis=2) which of course is a loop but it's done in the compiled code that numpy calls.

There may be some speed optimization possible here by adjusting the array definitions so that we can sum over a different axis and other things as well. See Code Review answers to

I've added random weights for demonstration purposes.

Convolute a gaussian kernel with a large array of off-grid centroids without looping?

might have been a somewhat misleading title as the problem is truly a sum over a finite number of individual centroids. The linked (and currently unanswered) question will be more of a challenge since it is a true convolution. 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
Yb = Y[..., None] - P
sig1 = 5
A = np.exp(-(Xb**2 + Yb**2) * sig1)
weights = np.random.random(len(P))
img = (A * weights).sum(axis=2)


The classical way to do this fast for arbitrary collections of "source" and "target" points is to use a fast multipole-type algorithm called the Fast Gauss Transform, developed by Greengard. A quick Google search turned up this Julia package implementing what the author deems an "improved" FGT: https://github.com/jwmerrill/FastGaussTransforms.jl

Based on "Improved Fast Gauss Transform," Proceedings International Conference on Computer Vision, 464 -471, C. Yang, R. Duraiswami, N.A. Gumerov, L. Davis. pdf

• I'll take a look, thanks!
– uhoh
Aug 16 at 13:59
• Ah, this is exactly what I want! It will take me some time to understand how to implement it in Python but I think it will be worth it. If I fail, I'll ask about it in a new question. Thanks!
– uhoh
Aug 16 at 15:38
• This appears to implement the FGT in Python: github.com/fredRos/pyfigtree Aug 16 at 19:13