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

user_14046024 0
Viewed 551

I need to draw a circle in a 2D numpy array given [i,j] as indexes of the array, and r as the radius of the circle. Each time a condition is met at index [i,j], a circle should be drawn with that as the center point, increasing all values inside the circle by +1. I want to avoid the for-loops at the end where I draw the circle (where I use p,q to index) because I have to draw possibly millions of circles. Is there a way without for loops? I also don't want to import another library for just a single task.

Here is my current implementation:

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
numpy geometry arrays python
2 Answers

If you're using the same radius for every single circle, you can simplify things significantly by only calculating the circle coordinates once and then adding the center coordinates to the circle points when needed. Here's the code:

# 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

When I set radius = 10, array = np.random.rand(1200).reshape(40, 30) and replaced if CONDITION with if (x == 20 and y == 20) or (x == 25 and y == 20), I got this, which seems to be what you want: Overlapping circles Let me know if you have any questions.

user_172292140

Adding each circle can be vectorized. This solution iterates over the coordinates where the condition is met. On a 2-core colab instance ~60k circles with radius 30 can be added per second.

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]

Visualizing np.array arr

import matplotlib.pyplot as plt

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

result

user_142777220
Related