How can I generate random points faster in a large arrays?

I am trying to generate random points inside a 3d pepper like shape however when the arr_size is large it takes too long to generate these points.

arr_size is (30,30,30) it takes very little time to generate 1000 random point inside the 3d shape however when arr_size = (265,490,286) it takes a realy long time.

from matplotlib import pyplot as plt
import numpy as np


def create_bin_pepper(arr_size, center):
    coords = np.ogrid[:arr_size[0], :arr_size[1], :arr_size[2]]
    c = 10
    a1 = np.random.randint(low=5,high=10)
    b1 = np.random.randint(low=7,high=10)
    a2 = np.random.randint(low=5,high=10)
    b2 = np.random.randint(low=7,high=10)
    a3 = np.random.randint(low=5,high=10)
    b3 = np.random.randint(low=7,high=10)
    ellipse1 = ((np.square(coords[0] - center[0]))/np.square(a1) + (np.square(coords[1]-center[1]))/np.square(b1)  + (np.square(coords[2]-center[2]))/np.square(c)  <= 1)
    ellipse2 = ((np.square(coords[0] - center[0]-3))/np.square(a2) + (np.square(coords[1]-center[1]-5))/np.square(b2)  + (np.square(coords[2]-center[2]))/np.square(c)  <= 1)
    ellipse3 = ((np.square(coords[0] - center[0]+3))/np.square(a3) + (np.square(coords[1]-center[1]-5))/np.square(b3)  + (np.square(coords[2]-center[2]))/np.square(c)  <= 1)
    pepper = ellipse1|ellipse2|ellipse3
    pepper2 = np.where(pepper==1,230,pepper)

    for im in range(0,1000):
        #r2=1
        centre_x1 = np.random.randint(low=center[0]-a1+4,high=center[0]+a1-4)#low=11,high=20
        centre_y1 = np.random.randint(low=center[1]-b1+4,high=center[1]+b1-4)#low=15,high=23
        centre_z1 = np.random.randint(low=center[2]-c+4,high=center[2]+c-4)#low=10,high=20

        centre_x2 = np.random.randint(low=center[0]-a2+4,high=center[0]+a2-4)#low=11,high=20
        centre_y2 = np.random.randint(low=center[1]-b2+4,high=center[1]+b2-4)#low=15,high=23
        centre_z2 = np.random.randint(low=center[2]-c+4,high=center[2]+c-4)#low=10,high=20

        centre_x3 = np.random.randint(low=center[0]-a3+4,high=center[0]+a3-4)#low=11,high=20
        centre_y3 = np.random.randint(low=center[1]-b3+4,high=center[1]+b3-4)#low=15,high=23
        centre_z3 = np.random.randint(low=center[2]-c+4,high=center[2]+c-4)#low=10,high=20

        inside_ellipse1 = ((np.square(coords[0] - centre_x1))/np.square(a1) + (np.square(coords[1]-centre_y1))/np.square(b1)  + (np.square(coords[2]-centre_z1))/np.square(c)  <= (1/((np.square(a1))*(np.square(b1))*(np.square(c)))))
        inside_ellipse2 = ((np.square(coords[0] - centre_x2-3))/np.square(a2) + (np.square(coords[1]-centre_y2-5))/np.square(b2)  + (np.square(coords[2]-centre_z2))/np.square(c)  <= (1/((np.square(a2))*(np.square(b2))*(np.square(c)))))
        inside_ellipse3 = ((np.square(coords[0] - centre_x3+3))/np.square(a3) + (np.square(coords[1]-centre_y3-5))/np.square(b3)  + (np.square(coords[2]-centre_z3))/np.square(c)  <= (1/((np.square(a3))*(np.square(b3))*(np.square(c)))))

    pepper2 = inside_ellipse1 | inside_ellipse2 | inside_ellipse3 | pepper2
    pepper3 = np.where((pepper2!=230)&(pepper2!=0),160,pepper2)
    return pepper3

arr_size = (265,490,286)
sphere_center1 = (133,216,40)

pepper = create_bin_pepper(arr_size,sphere_center1)
axis = pepper[:,:,40]
plt.imshow(axis,cmap='gray')#,interpolation='bicubic'
plt.show()