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.

"A Thousand (Gaussian) Points of Light"

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)

3 Answers 3


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

  • $\begingroup$ 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! $\endgroup$
    – uhoh
    Aug 16, 2021 at 15:38
  • 1
    $\begingroup$ This appears to implement the FGT in Python: github.com/fredRos/pyfigtree $\endgroup$ Aug 16, 2021 at 19:13
  • 1
    $\begingroup$ I think I am finally about to get my round tuit and learn/understand/use FGT, yay! :-) $\endgroup$
    – uhoh
    Mar 11, 2023 at 11:40

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
  • $\begingroup$ 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... $\endgroup$
    – uhoh
    Jul 21, 2021 at 13:26
  • 1
    $\begingroup$ 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)$ $\endgroup$
    – Ron
    Jul 21, 2021 at 13:29
  • $\begingroup$ Thanks, I've added an answer for Python based on your answer. $\endgroup$
    – uhoh
    Aug 15, 2021 at 4:39
  • $\begingroup$ I'm still puzzling over the mention of the discretization of a Fredholm integral equation of the first kind. Is this offered only as some helpful insight, or is this an avenue for a (potentially) substantially faster solution if there's some level of acceptable error, say ppt or ppm. (cf. Computer graphics III – Rendering equation and its solution) I think there can be a new question here (but first I'll check to see if this has already been covered) $\endgroup$
    – uhoh
    Dec 15, 2021 at 21:34

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.

enter image description here

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) 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.