python - 在给定索引和半径的 numpy 数组中绘制一个圆，无需外部库

Question

我需要在 2D numpy 数组中画一个圆，给定 [i,j] 作为数组的索引，r 作为圆的半径。每次在索引 [i,j] 处满足条件时，应以该点为中心点绘制一个圆，并将圆内的所有值增加 +1。我想避免在画圆的最后出现 for 循环(我使用 p,q 来索引)，因为我必须画可能数百万个圆。有没有不用for循环的方法？我也不想仅为一个任务导入另一个库。

这是我当前的实现:

for i in range(array_shape[0]):
    for j in range(array_shape[1]):
        if (condition):  # Draw circle if condition is fulfilled

            # Create a square of pixels with side lengths equal to radius of circle
            x_square_min = i-r
            x_square_max = i+r+1
            y_square_min = j-r
            y_square_max = j+r+1

            # Clamp this square to the edges of the array so circles near edges don't wrap around
            if x_square_min < 0:
                x_square_min = 0
            if y_square_min < 0:
                y_square_min = 0
            if x_square_max > array_shape[0]:
                x_square_max = array_shape[0]
            if y_square_max > array_shape[1]:
                y_square_max = array_shape[1]

            # Now loop over the box and draw circle inside of it
            for p in range(x_square_min , x_square_max):
                for q in range(y_square_min , y_square_max):
                    if (p - i) ** 2 + (q - j) ** 2 <= r ** 2:
                        new_array[p,q] += 1  # Incrementing because need to have possibility of 
                                             # overlapping circles              

Answer 1

添加每个圆都可以矢量化。此解决方案迭代满足条件的坐标。上2-core colab instance每秒可添加约 60k 个半径为 30 的圆。

import numpy as np

np.random.seed(42)
arr = np.random.rand(400,300)
r = 30

xx, yy = np.mgrid[-r:r+1, -r:r+1]
circle = xx**2 + yy**2 <= r**2

condition = np.where(arr > .999) # np.where(arr > .5) to benchmark 60k circles
for x,y in zip(*condition):
    # valid indices of the array
    i = slice(max(x-r,0), min(x+r+1, arr.shape[0]))
    j = slice(max(y-r,0), min(y+r+1, arr.shape[1]))

    # visible slice of the circle
    ci = slice(abs(min(x-r, 0)), circle.shape[0] - abs(min(arr.shape[0]-(x+r+1), 0)))
    cj = slice(abs(min(y-r, 0)), circle.shape[1] - abs(min(arr.shape[1]-(y+r+1), 0)))
    
    arr[i, j] += circle[ci, cj]

可视化np.array arr

import matplotlib.pyplot as plt

plt.figure(figsize=(8,8))
plt.imshow(arr)
plt.show()