Simulating Particle Image

simulate.py

def core_shell_gaussian(x, core, shell):
    """
    since we're dealing with REAL image with normally uint8 float type
    value < 1/255 --> value = 0
    let zero = 1/256
    """
    constant = shell ** 2 * 0.18033688011112042  # - shell^2 / ln(zero)
    core_shell = np.exp(- (x - core) ** 2 / constant)
    core_shell[x <= core] = 1.0  # construct core
    core_shell[x >= (core + shell)] = 0.0  # trim 4 edges of the box
    if core_shell.any() < 0:
        raise ValueError('wrong core-gaussian particle: negative intensity')
    return core_shell


def py_simulate_nd(image, positions, core, shell):
    new_image = np.zeros(image.shape, dtype=np.float64)
    dimensions, particle_numbers = positions.shape
    lower_bounds = np.zeros(positions.shape, dtype=int)
    upper_bounds = np.ones(positions.shape) * np.array([image.shape], dtype=int).T
    box_size = np.array([core + shell + 1] * 3)

    for dim in range(dimensions):
        lower_bounds[dim] = np.max([
            lower_bounds[dim],
            np.floor(positions[dim] - box_size[dim])
        ], axis=0)

    for num in range(particle_numbers):
        sub_image_box, sub_image_indices = [], []
        for dim in range(dimensions):
            sub_image_indices.append(
                np.linspace(
                    lower_bounds[dim, num] - positions[dim, num],
                    upper_bounds[dim, num] - positions[dim, num],
                    upper_bounds[dim, num] - lower_bounds[dim, num] + 0,
                    endpoint=True
                )
            )
            sub_image_box.append(
                slice(int(lower_bounds[dim, num]), int(upper_bounds[dim, num] + 1))
            )
        mesh = np.meshgrid(*sub_image_indices, indexing='ij', sparse=True)
        radius = np.sqrt(np.sum(np.array(mesh) ** 2, axis=0))
        if 0 in radius.shape:
            raise ValueError(f"Box for particle #{num} is outside the image")
        new_image[tuple(sub_image_box)] += core_shell_gaussian(radius, core, shell)
    return new_image