Размещение точек с переменной плотностью

Предположим, что у вас есть матрица NxM со значениями в диапазоне от [0,100]. Я бы хотел разместить точки с плотностью (обратно) относительно значений в этой области.

Например, вот двумерное гауссово поле, инвертированное s.t. центроид имеет значение 0, а периметр равен 100:

Я хотел бы упаковать точки так, чтобы они выглядели как на этом изображении:

Обратите внимание на радиальное распространение наружу.

Моя попытка выглядит немного иначе :( ...

Я пытаюсь (i) создать логическую область той же формы и размера и (ii) перемещаться по строкам и столбцам. Если значение логического массива в какой-то момент равно True, то передать; в противном случае добавьте точку [row,col] в список и закройте логический массив True в радиусе, пропорциональном значению в гауссовском массиве.

Выбор гауссиана для этого примера не важен, основная идея заключается в следующем: учитывая матрицу с плавающей запятой, как можно разместить точки с плотностью, пропорциональной этим значениям?

Любая помощь очень ценится :)

import matplotlib.pyplot as plt
import numpy as np
from math import exp

def gaussian(x,y,x0,y0,A=10.0,sigma_x=10.0,sigma_y=10.0):
    return A - A*exp(-((x-x0)**2/(2*sigma_x**2) + (y-y0)**2/(2*sigma_y**2)))

def generate_grid(width=100,height=100):
    grid = np.empty((width,height))
    for x in range(0,width):
        for y in range(0,height):
            grid[x][y] = gaussian(x,y,width/2,height/2,A=100.0)
    return grid

def cover_array(a,row,col,radius):
    nRows = np.shape(grid)[0]
    nCols = np.shape(grid)[1]
    mid = round(radius / 2)
    half_radius = int(round(radius))
    for x in range(-half_radius,half_radius):
        for y in range(-half_radius,half_radius):
            if row+x >= 0 and x+row < nRows and col+y >= 0 and y+col < nCols:
                if (x-mid)**2 + (y-mid)**2 <= radius**2:
                    a[row+x][col+y] = True


def pack_points(grid):
    points = []
    nRows = np.shape(grid)[0]
    nCols = np.shape(grid)[1]
    maxDist = 50.0
    minDist = 0.0
    maxEdge = 10.0
    minEdge = 5.0
    grid_min = 0.0
    grid_max = 100.0

    row = 0
    col = 0

    arrayCovered = np.zeros((nRows,nCols))

    while True:
        if row >= nRows:
            return np.array(points)

        if arrayCovered[row][col] == False:
            radius = maxEdge * ((grid[row][col] - grid_min) / (grid_max - grid_min))
            cover_array(arrayCovered,row,col,radius)
            points.append([row,col])

        col += 1
        if col >= nCols:
            row += 1
            col = 0

grid = generate_grid()
plt.imshow(grid)
plt.show()

points = pack_points(grid)

plt.scatter(points[:,0],points[:,1])
plt.show()