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
- Better way to calculate double-scattering diffraction using cartesian product of arrays?
- improving speed of this numpy-based diffraction calculator
- Getting hexagonal arrays of dots arranged in this spiral pattern
- and particularly Calculate electric field of a charged annulus
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[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)