My attempt would look something like this:

if we're concerned about coverage, we can change it thusly:

But I don't have a good enough intuition for barycentric coordinates to know which of those we prefer without drawing it out better.

**Code:**```
```

indices=np.arange(3)

oboe = (indices+1)%3

# First get surface normal and edge lengths! :)

# e = tv[oboe] - tv[indices] # <- cool way to do it

e = array([tv[oboe[i]] - tv[i] for i in indices]) # <- boring way

#L1=numpy.linalg.norm(e1)

#L2=numpy.linalg.norm(e2)

#L3=numpy.linalg.norm(e3)

L = np.apply_along_axis(np.linalg.norm,-1,e)

snorm = numpy.cross(e[1],e[0])

snorm = snorm/numpy.linalg.norm(snorm)

Linc = 1 /L

Lspan = np.ceil(L)

Lspaces = [np.linspace(0,1,1+Lspan[i]) for i in indices]

bary_coords = np.zeros(3)

for i,j in numpy.vstack((indices,oboe)).transpose():

for ca in Lspaces[i]:

for cb in Lspace[j]:

if ca+cb > 1:

break

bary_coords.fill(1 - ca - cb)

bary_coords[i] = ca

bary_coords[j] = cb

# here we will need to have calculated matrices

# to take us from barycentric coords to 3-space

# Cast ray must take two args, src and dir

# and moves from src to src+dir

castraythroughvoxels(make_cartesian_coords(bary_coords),snorm)

if we're concerned about coverage, we can change it thusly:

**Code:**```
```

indices=np.arange(3)

oboe = (indices+1)%3

# First get surface normal and edge lengths! :)

# e = tv[oboe] - tv[indices] # <- cool way to do it

e = array([tv[oboe[i]] - tv[i] for i in indices]) # <- boring way

#L1=numpy.linalg.norm(e1)

#L2=numpy.linalg.norm(e2)

#L3=numpy.linalg.norm(e3)

L = np.apply_along_axis(np.linalg.norm,-1,e)

snorm = numpy.cross(e[1],e[0])

snorm = snorm/numpy.linalg.norm(snorm)

Linc = 1 /L

Lspan = np.ceil(L)

Lspaces = [np.linspace(0,1,1+Lspan[i]) for i in indices]

bary_coords = np.zeros(3)

for i,j in numpy.vstack((indices,oboe)).transpose():

for ca in Lspaces[i]:

for cb in Lspace[j]:

if ca+cb > 1:

break

bary_coords.fill(1 - ca - cb)

bary_coords[i] = ca

bary_coords[j] = cb

# here we will need to have calculated matrices

# to take us from barycentric coords to 3-space

# Cast ray must take two args, src and dir

# and moves from src to src+dir

castraythroughvoxels(make_cartesian_coords(bary_coords),snorm)

for ca in Lspaces[j][len(Lspaces[j])/2:]:

for cb in Lspace[i]:

if ca+cb > 1:

break

bary_coords.fill(1 - ca - cb)

bary_coords[j] = ca

bary_coords[i] = cb

# here we will need to have calculated matrices

# to take us from barycentric coords to 3-space

# Cast ray must take two args, src and dir

# and moves from src to src+dir

castraythroughvoxels(make_cartesian_coords(bary_coords),snorm)

But I don't have a good enough intuition for barycentric coordinates to know which of those we prefer without drawing it out better.