How to fill 3D figure with voxels?

I'm attempting to fill a tetrahedron with voxels to generate 3D data. I have been able to generate the tetrahedron itself using 4 different points. I'm not exactly sure how I could use NumPy or any other Python framework for filling the region inside the tetrahedron with voxels. Here is the code for generating the tetrahedron 3D plot:

# Tetrahedron

from scipy.spatial import Delaunay
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import numpy as np
points= np.array([[1,2,2],[1,3,4],[4,1,4],[5,3,2]])
tri = Delaunay(points)
tr = tri.simplices[0]  # indices of first tetrahedron
pts = points[tr, :]  # pts is a 4x3 array of the tetrahedron coordinates
fig = plt.figure()
ax = fig.gca(projection= '3d')

# plotting the six edges of the tetrahedron
for ij in [[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]]:
    ax.plot3D(pts[ij, 0], pts[ij, 1], pts[ij, 2])
plt.show()

Here is an example image of what I am trying to achieve. This example features a sphere that has been filled with voxels:

So you're looking for something like this?

I have adopted this https://matplotlib.org/stable/gallery/mplot3d/voxels_rgb.html to get your simplex.

import itertools
import functools
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["figure.figsize"] = (18,9)
points = np.array([[1,2,2],[1,3,4],[4,1,4],[5,3,2]])
center = np.mean(points, axis=0)

def surface_normal_form(a,b,c):
    v = b-a
    w = c-b
    n = np.cross(v,w)
    #normal needs to point out
    if (center-a)@n > 0:
         n *= -1
    return a, n

def midpoints(x):
    sl = ()
    for i in range(x.ndim):
        x = (x[sl + np.index_exp[:-1]] + x[sl + np.index_exp[1:]]) / 2.0
        sl += np.index_exp[:]
    return x

x, y, z = (np.indices((60, 60, 60))-np.array([20,25,25]).reshape(-1,1,1,1))/8
mx = midpoints(x)
my = midpoints(y)
mz = midpoints(z)

conditions = []
for p1,p2,p3 in itertools.combinations(points, 3):
    a, n = surface_normal_form(p1,p2,p3)
    conditions.append((mx-a[0])*n[0]+(my-a[1])*n[1]+(mz-a[2])*n[2] <= 0)

simplex = conditions[0] & conditions[1] & conditions[2] & conditions[3]

ax = plt.figure().add_subplot(projection='3d')
ax.voxels(x, y, z, simplex, linewidth=0.5)
ax.set(xlabel='x', ylabel='y', zlabel='z')
ax.set_xlim(1.0,5.0)
ax.set_ylim(1.0,3.0)
ax.set_zlim(2.0,4.0)

Can you explain what you're doing in the surface_normal_form function?

Sure in 3d you can describe a plane by a point of the plane and vector that's orthogonal to the plane. That's particularly useful because it's then easy to tell if a point is on one side of the plane or the other. If you then take the planes containing the sides of the simplex you get get the simplex voxels by taking a cube of voxels tightly fit together and for each plane removing the ones on the wrong side of it. Thats what I am doing.

Also, can you please explain the midpoints function as well?

First a little disclaimer I haven't written that myself. As I said it's from the matplotlib example. But if you wanted to calculate midpoints of a 1d array you could do.

(x[:-1]+x[1:])/2

Firstly x[:-1] gives you all the values but the last and x[1:] gives you all the values but the first so the values that are next to each other getting summed and divided by two i.e. you get the midpoints. Notice that it takes the "midpoints" of (in our case 2d) subarrays if the original array is more than 1d (in our case it's 3 dimensional). In step two what we're doing is (x[:,:-1]+x[:,1:])/2. And since [:] gives you all the values the midpoint function does that for each dimension.