A lot of calculations can be simplified by turning this into a 2D problem, and you can do that by casting your ray along a single axis -- x, y, or z.

For a given V = { Vx, Vy, Vz }, if we cast our ray along the z axis, we're effectively just checking which triangles intersect the 2D point V2 = { Vx, Vy } (projected into 2D by removing the z component). From there, find the z value for the point of intersection of V with each triangle. Count the number of triangles whose z-point of intersection is greater than Vz (we're effectively doing a ray trace in the positive z direction). If it's even, V is outside the mesh; if it's odd, it's inside the mesh (assuming the mesh is closed).

This leads to the main optimization: Rasterise the triangles into a 2D grid.

That is, you'll want to pre-calculate which triangles intersect which 2D { x, y } points in your grid, and with what z value, because that way you're only calculating triangle intersections for 128^2 points instead of 128^3 points. This can be done by rasterising the triangles to a 2D 128^2 grid, but keeping a list of potentially multiple values for each point in the grid. The 2D grid could be anything from a 2D array of arrays to something more space-efficient, but I think that's beyond the scope of this question.

A lot of calculations can be simplified by turning this into a 2D problem, and you can do that by casting your ray along a single axis -- x, y, or z.

For a given V = { Vx, Vy, Vz }, if we cast our ray along the z axis, we're effectively just checking which triangles intersect the 2D point V2 = { Vx, Vy } (projected into 2D by removing the z component). From there, find the z value for the point of intersection of V with each triangle. Count the number of triangles whose z-point of intersection is greater than Vz (we're effectively doing a ray trace in the positive z direction). If it's even, V is outside the mesh; if it's odd, it's inside the mesh (assuming the mesh is closed).

This leads to the main optimization: Rasterise the triangles into a 2D grid.

That is, you'll want to pre-calculate which triangles intersect which 2D { x, y } points in your grid, and with what z value, because that way you're only calculating triangle intersections for 128^2 points instead of 128^3 points. This can be done by rasterising the triangles to a 2D 128^2 grid, but keeping a list of intersecting triangles or z points of intersection for each point in the grid. The 2D grid could be anything from a 2D array of arrays to something more space-efficient, but I think that's beyond the scope of this question.