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) 

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.

• How to generate the convolution of f(x, y) with a parametric function g(t), x(t), y(t) in Python? (Something better than this brute-force sum)

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)