How to fill selected cells in 3D Matplotlib?

This should a straight forward question but I could not find the answer, sorry if the question is a duplicate one. Basically I want to fill Ncell cells on a 3D grid of Ngrid x Ngrid x Ngrid. Below I present a MWE where I just scatter plot the centre (or one corner, it does not matter, I can readjust) of each cell.

import numpy as np
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt

Ngrid,Ncell=100,150

xx=np.random.choice(Ngrid,Ncell)
yy=np.random.choice(Ngrid,Ncell)
zz=np.random.choice(Ngrid,Ncell)

%matplotlib widget

fig = plt.figure(figsize = (10, 7))
ax = plt.axes(projection ="3d")
ax.scatter3D(xx,yy,zz,marker='o',alpha=0.5)

Instead of the scatter plot, I need each cell to be filled with a color. Note that my grid and Ncell are much larger than the above values (grid = 1000^3 and Ncell= order 10^5), so an efficient code will be very important. Thank you in advacne.

You can try with voxels but it is already slow with your test data:

import numpy as np
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt

Ngrid,Ncell=100,250

xx=np.random.choice(Ngrid,Ncell)
yy=np.random.choice(Ngrid,Ncell)
zz=np.random.choice(Ngrid,Ncell)
v = np.zeros((Ngrid, Ngrid, Ngrid), dtype=bool)
v[xx, yy, zz] = True

fig = plt.figure(figsize = (10, 7))
ax = plt.axes(projection ="3d")
ax.voxels(v)

I would suggest to use K3D Jupyter and its voxel. This appears to be much much faster. Note that, as the name suggests, it requires Jupyter Notebook.

import k3d
import numpy as np

Ngrid,Ncell=100,1000

xx=np.random.choice(Ngrid,Ncell)
yy=np.random.choice(Ngrid,Ncell)
zz=np.random.choice(Ngrid,Ncell)

v = np.zeros((Ngrid, Ngrid, Ngrid), dtype=bool)

# NOTE: K3D's voxels works with the the matrix indexing, ij, rather
# then with the cartesian indexing, xy. Read np.meshgrid and run
# some example to understand the differences.
# The bottom line is, the following is the correct indexing:
v[zz, yy, xx] = True

# Once the plot is shown on the screen we can play with options
# by clicking K3D panel -> Objects -> Voxels #1. If we find some
# good options, we can write the values in this command and
# execute it again.
plt_voxels = k3d.voxels(
    v.astype(np.uint8),
    color_map=(0xfdfe03),
    outlines=False, # better hide outlines when setting opacity < 1
    opacity=0.5
)
plot = k3d.plot(grid_visible=False)
plot += plt_voxels

# create scatter
positions = np.vstack([t.flatten() for t in [xx+0.5, yy+0.5, zz+0.5]]).T.astype(np.float32)
plt_points = k3d.points(positions=positions, point_size=0.4, color=0xff0000)
plot += plt_points

plot.display()