That function looks pretty well behaved/smooth. So assuming your voxels are pretty small compared to the surface curvature, I think you could get away with a signed-distance-like approach, by testing the 8 corners against the surface. This will quickly eliminate the all-corners-inside (1) and all-corners-outside (0) cases.
Near the surface there will be ...
If you are familiar with einsum, maybe this explanation does it: axes and axes specify the locations of the repeated letters in the parameters of einsum. For instance,
np.tensordot(a, b, axes=[(0,2),(3,1)])
np.einsum('ijkl,mkni', a, b)
Indeed, 'ijkl'[(0,2)] == 'ik' == 'mkni'[(3,1)], and all the other letters are distinct.