Fill a numpy array with ones within a circle, zeros elsewhere

I am trying to imitate a circle filled with ones and the rest of the array filled with zeros in a numpy array. This has something to do with a numerical computation of the overlap of a figure with circles. For now at best I can do it elementwise (which is slow):

import numpy as np
from my_magic_constants import center, radius

c1 = np.zeros((center*2, center*2))
for h in range(center*2):
	for w in range(center*2):
		if (center-h)**2 + (center-w)**2 < radius**2:
			c1[h,w] = 1

What is the more efficient way to initialize a circle in a numpy array?

Update

I have only checked the acceleration by the earliest answer. If someone wants to benchmark timings, here is the code

import numpy as np
from time import time

eps = 1e-10
center = 8119
square_half_size = 5741
square = np.zeros((center*2, center*2))
square[center-square_half_size:center+square_half_size, center-square_half_size:center+square_half_size] = 1


def iou(c):
    overlap = c * square
    union = (c + square).clip(0,1)
    return np.sum(overlap)/(np.sum(union) + eps)


best_iou = 0

for radius in range(center - int(square_half_size * 1.093), center):
    t1 = time()
    c1 = np.zeros((center*2, center*2))
    for h in range(center*2):
        for w in range(center*2):
            if (center-h)**2 + (center-w)**2 < radius**2:
                c1[h,w] = 1

    current_iou = iou(c1)
    if current_iou > best_iou:
        best_iou = current_iou
        # print(i, best_iou)
    t = time() - t1

    print(t)
    break

If some later answer is 10+ times faster, we should acknowledge that.

You can also use broadcasting, or somewhat equivalently np.add.outer:

import numpy as np
center, radius = 10, 5

arr = np.arange(-center, center) ** 2
out = np.add.outer(arr, arr) < radius ** 2
# or: arr[:, None] + arr[None, :] < radius ** 2

# Optionally convert dtype if you don't like bool:
out = out.astype(float)

You can use a Boolean mask. The key is to have a way to create all h and w at the same time, below called X and Y. Check out the meshgrid docs for details.

import numpy as np
center, radius = 10, 5

c1 = np.zeros((center*2, center*2))
for h in range(center*2):
    for w in range(center*2):
        if (center-h)**2 + (center-w)**2 < radius**2:
            c1[h,w] = 1
print(c1)

c2 = np.zeros((center*2, center*2))
x = np.arange(c2.shape[0])
X, Y = np.meshgrid(x, x)
c2[(X - center)**2 + (Y - center)**2 < radius**2] = 1
print(c2)

print((c1 == c2).all())  # True