0

The convolution of a regular grid of Dirac delta functions with a kernel is pretty standard:

import numpy as np
deltas = np.zeros((2048, 2048))
deltas[8::16,8::16] = 1

# Construct a Gaussian kernel
x, y = np.meshgrid(np.linspace(-2, 2, 15), np.linspace(-2, 2, 15))
gauss = np.exp(-(x*x + y*y)/2)
from scipy.signal import convolve2d
convolved_image = convolve2d(deltas, gauss)

Then you get all your Gaussians on a regular grid. An FFT convolution would be a lot faster and is as easy to implement, but that is not my issue.

Now, in order to generate input data for a statistical test, the Gaussians should not be centered exactly at the middle of a pixel, but have a random position anywhere on a pixel. So the chance of its peak being positioned in the lower left corner of a pixel should be equal to its peak being positioned at the center of pixel.

Besides these random positional fluctuations at the pixel level, I would still like to produce a regular grid of Gaussian, so from a distance it would still look like a perfectly regular grid.

How would you code this in Python?

Alex van Houten
  • 308
  • 3
  • 17
  • A Gaussian with a random uniform offset is no longer a Gaussian, the pdf is the convolution of a Gaussian with the uniform distribution. You can compute this analytically (difference of two error functions) and use that as your input. – codeMonkey Jun 19 '23 at 07:13
  • I would call f = a * exp( (x-b)**2/(2 *c**2)), which has an offset b, a Gaussian. According to Wikipedia this is a Gaussian https://en.wikipedia.org/wiki/Gaussian_function – Alex van Houten Jun 20 '23 at 18:47
  • I know the definition of a Gaussian ;) This is not the distribution (kernel) you seek if "the chance of its peak being positioned in the lower left corner of a pixel should be equal to its peak being positioned at the center of pixel". This distribution is flat topped. You need to convolve the Gaussian with a uniform distribution analytically to get that distribution if you do not want to do it numerically like in Serge's answer below. – codeMonkey Jun 27 '23 at 11:06
  • Ah, so what you mean is that there is an analytical solution to this problem? – Alex van Houten Jun 29 '23 at 08:44

1 Answers1

0

In order to do this you can add a random offset to each Gaussian function by shifting the grid before applying the convolution. Here the random_offset is calculated once and then added to both x and y grids. The same offset is used for every Gaussian function in the grid.

import numpy as np
import random

deltas = np.zeros((2048, 2048))
deltas[8::16,8::16] = 1

x, y = np.meshgrid(np.linspace(-2, 2, 15), np.linspace(-2, 2, 15))
gauss = np.exp(-(x*x + y*y)/2)

# add a random offset to the grid
random_offset = (random.uniform(-0.5, 0.5), random.uniform(-0.5, 0.5))
x = x + random_offset[0]
y = y + random_offset[1]
gauss = np.exp(-(x*x + y*y)/2)

from scipy.signal import convolve2d
convolved_image = convolve2d(deltas, gauss)
Serge de Gosson de Varennes
  • 7,162
  • 3
  • 18
  • 39
  • Thanks. Nice way to start, but this could become a long run, I reckon, since you would have to loop over (2048/16)² = 16384 convolutions. And then add them all up. Perhaps this cannot be avoided, but making deltas = np.zeros((16, 16)); deltas[8, 8] = 1, would make this somewhat faster. Would require aggregating the convolved subimages into a 2048 * 2048 final image. – Alex van Houten Feb 13 '23 at 10:28