2
$\begingroup$

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)
    plt.show()
$\endgroup$
1
4
$\begingroup$

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)
$\endgroup$
3

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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