在没有外部库的情况下,在给定索引和半径的 numpy 数组中绘制一个圆

Draw a circle in a numpy array given index and radius without external libraries

我需要在二维 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              

如果您对每个圆都使用相同的半径,则可以通过仅计算一次圆坐标然后在需要时将圆心坐标添加到圆点来显着简化事情。这是代码:

# The main array of values is called array.
shape = array.shape
row_indices = np.arange(0, shape[0], 1)
col_indices = np.arange(0, shape[1], 1)

# Returns xy coordinates for a circle with a given radius, centered at (0,0).
def points_in_circle(radius):
  a = np.arange(radius + 1)
  for x, y in zip(*np.where(a[:,np.newaxis]**2 + a**2 <= radius**2)):
    yield from set(((x, y), (x, -y), (-x, y), (-x, -y),))

# Set the radius value before running code.
radius = RADIUS
circle_r = np.array(list(points_in_circle(radius)))

# Note that I'm using x as the row number and y as the column number.
# Center of circle is at (x_center, y_center). shape_0 and shape_1 refer to the main array
# so we can get rid of coordinates outside the bounds of array.
def add_center_to_circle(circle_points, x_center, y_center, shape_0, shape_1):
  circle = np.copy(circle_points)
  circle[:, 0] += x_center
  circle[:, 1] += y_center

  # Get rid of rows where coordinates are below 0 (can't be indexed)
  bad_rows = np.array(np.where(circle < 0)).T[:, 0]
  circle = np.delete(circle, bad_rows, axis=0)

  # Get rid of rows that are outside the upper bounds of the array.
  circle = circle[circle[:, 0] < shape_0, :]
  circle = circle[circle[:, 1] < shape_1, :]

  return circle

for x in row_indices:
  for y in col_indices:

    # You need to set CONDITION before running the code.
    if CONDITION:
      # Because circle_r is the same for all circles, it doesn't need to be recalculated all the time. All you need to do is add x and y to circle_r each time CONDITION is met.
      circle_coords = add_center_to_circle(circle_r, x, y, shape[0], shape[1])

      array[tuple(circle_coords.T)] += 1

当我设置 radius = 10array = np.random.rand(1200).reshape(40, 30) 并将 if CONDITION 替换为 if (x == 20 and y == 20) or (x == 25 and y == 20) 时,我得到了这个,这似乎是你想要的: 如果您有任何问题,请告诉我。

添加每个圆都可以矢量化。此解决方案迭代满足 条件 的坐标。在 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()