Evaluating a general volume integral exactly, in Cartesian coordinates, on a grid of cubic voxels is relatively easy:
$$ G_i = \int_{X_{min}}^{X_{max}}\int_{Y_{min}}^{Y_{max}}\int_{Z_{min}}^{Z_{max}} g(x,y,z) dxdydz$$
where the limits are simply found from the max and mins of the set of coordinates defining each of the voxels.
How could this be approached for an arbitrary tetrahedral voxel instead, and how would the limits be chosen from the four vertex points assuming $g(x,y,z)$ is a known function?
(Use voxel, not pixel.)
If you discretize the volume into arbitrary-sized and -shaped cells (tetrahedra or whatever else), then the integral becomes a sum. For example, integrating a function $g$ over some specific volume discretized into cells, $$\iiint g(x,y,z) ~ d x ~ d y ~ d z = \iiint r^2 \sin(\theta) ~ g(r, \theta, \varphi) ~ d r ~ d \theta ~ d \varphi = \sum_i G_i$$ where $G_i$ is or approximates the volume integral over cell $i$. Typical approximations use $G_i = g_i V_i$ where $g_i$ is the average value of $g$ within the cell, and $V_i$ is the volume of the cell.
If cell $i$ is tetrahedral with vertices $\vec{a}_i$, $\vec{b}_i$, $\vec{c}_i$, and $\vec{d}_i$, then its volume $V_i$ is $$V_i = \frac{1}{6}\left\lvert(\vec{a}_i-\vec{d}_i)\cdot(\vec{b}_i-\vec{d}_i)\times(\vec{c}_i-\vec{d}_i)\right\rvert$$ where the factor $1/6$ is the factor OP wondered about. (Without the factor, $V_i$ would be the volume of a cube, where $\vec{d}_i$ is one of the eight vertices, and $\vec{a}$, $\vec{b}$, and $\vec{c}$ are the three neighboring vertices sharing an edge with vertex $\vec{d}$.
If I understood OP correctly, the problem is essentially how to calculate the exact integral for $G_i$ when $g(x,y,z)$ and the four vertices of the bounding tetrahedron are known, $$G_i = \iiint_{V_i} g(x, y, z) ~ d x ~ d y ~ d z$$
The tetrahedron is bounded by four planar faces. Each plane is defined by its normal vector $\vec{n}$ and signed distace from origin $d$: point $\vec{p}$ is on the plane if and only if $$\vec{n}\cdot\vec{p} = d$$ The normal vectors can be computed from the vertices trivially: $$\begin{aligned} \vec{n}_{abc} &= (\vec{b} - \vec{a})\times(\vec{c} - \vec{a}) \\ \vec{n}_{abd} &= (\vec{b} - \vec{a})\times(\vec{d} - \vec{a}) \\ \vec{n}_{acd} &= (\vec{c} - \vec{a})\times(\vec{d} - \vec{a}) \\ \vec{n}_{bcd} &= (\vec{c} - \vec{b})\times(\vec{d} - \vec{b}) \\ \end{aligned}$$ and can the signed distances: $$\begin{aligned} d_{abc} &= \vec{n}_{abc} \cdot \vec{a} = \vec{n}_{abc} \cdot \vec{b} = \vec{n}_{abc} \cdot \vec{c} \\ d_{abd} &= \vec{n}_{abd} \cdot \vec{a} = \vec{n}_{abd} \cdot \vec{b} = \vec{n}_{abd} \cdot \vec{d} \\ d_{acd} &= \vec{n}_{acd} \cdot \vec{a} = \vec{n}_{acd} \cdot \vec{d} = \vec{n}_{acd} \cdot \vec{c} \\ d_{bcd} &= \vec{n}_{bcd} \cdot \vec{b} = \vec{n}_{bcd} \cdot \vec{c} = \vec{n}_{bcd} \cdot \vec{d} \\ \end{aligned}$$ Note that each normal and the associated signed distance from origin can be multiplied by any nonzero scalar $\lambda$ without affecting the plane itself. Depending on the order of the vertices $\vec{a}$, $\vec{b}$, $\vec{c}$, and $\vec{d}$, the normal vectors can point either inwards or outwards; you can change that by multiplying both normal vector and signed distance by $-1$.
If you need the normals to point outwards, calculate centroid $\vec{p} = \frac{\vec{a} + \vec{b} + \vec{c} + \vec{d}}{4}$, and verify that $$\begin{aligned} \vec{n}_{abc} \cdot \vec{p} & \lt d_{abc} \\ \vec{n}_{abd} \cdot \vec{p} & \lt d_{abd} \\ \vec{n}_{acd} \cdot \vec{p} & \lt d_{acd} \\ \vec{n}_{bcd} \cdot \vec{p} & \lt d_{bcd} \\ \end{aligned}$$ (If not, negate both $\vec{n}$ and $d$ for that case.)
Using $$\vec{a} = \left [ \begin{matrix} a_x \\ a_y \\ a_z \end{matrix} \right ] , \quad \vec{b} = \left [ \begin{matrix} b_x \\ b_y \\ b_z \end{matrix} \right ] , \quad \vec{c} = \left [ \begin{matrix} c_x \\ c_y \\ c_z \end{matrix} \right ] , \quad \vec{d} = \left [ \begin{matrix} d_x \\ d_y \\ d_z \end{matrix} \right ]$$ we can write the integral in Cartesian coordinates is $$G = \int_{z_{min}}^{z_{max}} \int_{y_{min}(z)}^{y_{max}(z)} \int_{x_{min}(y, z}^{x_{max}(y, z)} g(x, y, z) ~ d x ~ d y ~ d z$$ where $z_{min} = \min(a_z, b_z, c_z, d_z)$, $z_{max} = \max(a_z, b_z, c_z, d_z)$; $y_{min}(z)$ and $y_{max}(z)$ give the smallest and largest $y$ coordinates within the tetrahedron at the given $z$ coordinate; and $x_{min}(y, z)$ and $x_{max}(y, z)$ give the smallest and largest $x$ coordinates within the tetrahedron. The problem is that the exact algebraic expressions for $y_{min}(z)$, $y_{max}(z)$, $x_{min}(y, z)$, and $x_{max}(y, z)$ are rather complex, unless we reorder the vertices.
Let's say we reorder the vertices according to their $z$ coordinates, so that $d_z \le c_z \le b_z \le a_z$. This also means that $z_{min} = d_z$ and $z_{max} = z_{max}$. Then, interpolating along each edge of the tetrahedron, $$\begin{aligned} \vec{e}_{da}(z) &= \vec{d} + \frac{z - d_z}{a_z - d_z} ( \vec{a} - \vec{d} ), \quad d_z \le z \le a_z \\ \vec{e}_{db}(z) &= \vec{d} + \frac{z - d_z}{b_z - d_z} ( \vec{b} - \vec{d} ), \quad d_z \le z \le b_z \\ \vec{e}_{dc}(z) &= \vec{d} + \frac{z - d_z}{c_z - d_z} ( \vec{c} - \vec{d} ), \quad d_z \le z \le c_z \\ \vec{e}_{ab}(z) &= \vec{a} + \frac{z - a_z}{b_z - a_z} ( \vec{b} - \vec{a} ), \quad b_z \le z \le a_z \\ \vec{e}_{ac}(z) &= \vec{a} + \frac{z - a_z}{c_z - a_z} ( \vec{c} - \vec{a} ), \quad c_z \le z \le a_z \\ \vec{e}_{bc}(z) &= \vec{b} + \frac{z - b_z}{c_z - b_z} ( \vec{c} - \vec{b} ), \quad c_z \le z \le b_z \\ \end{aligned}$$ At any given $z$, $d_z \le z \le a_z$, these yield the vertices of the triangle or quadrangle that is formed from the intersection of the plane $z$ and the tetrahedron. There are three regions: $$\begin{array}{c|l} z & \text{Valid edges} \\ \hline d_z \le z \le c_z & \vec{e}_{da}, \vec{e}_{db}, \vec{e}_{dc} \\ c_z \lt z \lt b_z & \vec{e}_{da}, \vec{e}_{db}, \vec{e}_{ac}, \vec{e}_{bc} \\ b_z \le z \le a_z & \vec{e}_{da}, \vec{e}_{ab}, \vec{e}_{ac} \\ \end{array}$$ so technically, even with the help of rearranging the vertices, $$y_{min}(z) = \begin{cases} \min([\vec{e}_{da}(z)]_y, [\vec{e}_{db}(z)]_y, [\vec{e}_{dc}(z)]_y), & d_z \le z \le c_z \\ \min([\vec{e}_{da}(z)]_y, [\vec{e}_{db}(z)]_y, [\vec{e}_{ac}(z)]_y, [\vec{e}_{bc}(z)]_y), & c_z \lt z \lt b_z \\ \min([\vec{e}_{da}(z)]_y, [\vec{e}_{ab}(z)]_y, [\vec{e}_{ac}(z)]_y), & b_z \le z \le a_z \\ \end{cases}$$ and same for $y_{max}(z)$ (except using $\max$ instead of $\min$).
For $x_{min}(y, z)$ and $x_{max}(y, z)$, we pick the three or four valid edge points based on $z$. For example, if $c_z \le z \le b_z$, the edge points are $\vec{e}_{da}(z)$, $\vec{e}_{db}(z)$, $\vec{e}_{ac}(z)$, and $\vec{e}_{bc}(z)$. In the $z$ plane, the four edges of the intersection are formed by line segments $\vec{e}_{da}(z),\vec{e}_{ac}(z)$, $\vec{e}_{db}(z),\vec{e}_{bc}(z)$, $\vec{e}_{ac}(z),\vec{e}_{bc}(z)$, and $\vec{e}_{da}(z),\vec{e}_{db}(z)$. The line segment is valid for a particular $y$ when one endpoint $\ge y$ and the other $\lt y$; then, the point $\vec{q}$ corresponding to $y$ on that line segment $\vec{e}_1,\vec{e}_2$ is $$\vec{q}(y) = \vec{e}_1 + \frac{y - [\vec{e}_1]_y}{[\vec{e}_2]_y - [\vec{e}_1]_y} ( \vec{e}_2 - \vec{e}_1 )$$ although at this point, only the $x$ coordinates of $\vec{q}$ matter: $x_{min}(y,z)$ and $x_{max}(y, z)$ are simply the smallest and largest $x$ coordinate among the points corresponding to $y$ on the valid line segments between the valid edges for $z$.
Another option is to use barycentric coordinates $(u, v, w)$ for the points within the tetrahedron. The four limiting planes are $u = 0$, $v = 0$, $w = 0$, and $u + v + w = 1$, and the integral becomes $$\iiint_{V_i} g(x, y, z) ~ d x ~ d y ~ d z = \int_{0}^{1} \int_{0}^{1-u} \int_{0}^{1-u-v} g\bigr(x(u, v, w), y(u, v, w), z(u, v, w)\bigr) V_d(u, v, w) ~ d u ~ d v ~ d w$$ where the differential volume unit $V_d$ is $$V_d(u, v, w) = \det \left [ \begin{matrix} \displaystyle \frac{\partial x(u,v,w)}{\partial u} & \displaystyle \frac{\partial x(u,v,w)}{\partial v} & \displaystyle \frac{\partial x(u,v,w)}{\partial w} \\ \displaystyle \frac{\partial y(u,v,w)}{\partial u} & \displaystyle \frac{\partial y(u,v,w)}{\partial v} & \displaystyle \frac{\partial y(u,v,w)}{\partial w} \\ \displaystyle \frac{\partial z(u,v,w)}{\partial u} & \displaystyle \frac{\partial z(u,v,w)}{\partial v} & \displaystyle \frac{\partial z(u,v,w)}{\partial w} \\ \end{matrix} \right ]$$ The inverse mapping from Barycentric to Cartesian coordinates, when the vertices of the tetrahedron are at Barycentric coordinates $\vec{d} = (0,0,0)$, $\vec{a} = (1,0,0)$, $\vec{b} = (0,1,0)$, and $\vec{c} = (0,0,1)$, is nonlinear $$\left\lbrace ~ \begin{aligned} x(u, v, w) &= (1 - w) \bigr( (1 - v) \bigr( (1 - u) d_x + u a_x \bigr) + v b_x \bigr) + w c_x \\ y(u, v, w) &= (1 - w) \bigr( (1 - v) \bigr( (1 - u) d_y + u a_y \bigr) + v b_y \bigr) + w c_y \\ z(u, v, w) &= (1 - w) \bigr( (1 - v) \bigr( (1 - u) d_z + u a_z \bigr) + v b_z \bigr) + w c_z \\ \end{aligned} \right .$$ The mapping from Cartesian to Barycentric coordinates is linear, $$\left\lbrace ~ \begin{aligned} u(x,y,z) &= u_0 + x u_x + y u_y + z u_z \\ v(x,y,z) &= v_0 + x v_x + y v_y + z v_z \\ w(x,y,z) &= w_0 + x w_x + y w_y + z w_z \\ \end{aligned} \right .$$ or equivalently in vector notation, $$\left\lbrace ~ \begin{aligned} u(\vec{p}) &= u_0 + \hat{u} \cdot \vec{p} \\ v(\vec{p}) &= v_0 + \hat{v} \cdot \vec{p} \\ w(\vec{p}) &= w_0 + \hat{w} \cdot \vec{p} \\ \end{aligned} \right .$$ where $$\begin{aligned} T &= \left(\vec{a}-\vec{d}\right) \cdot \bigr( \left(\vec{b} - \vec{d}\right) \times \left(\vec{c} - \vec{d}\right) \bigr) \\ u_0 &= \frac{ \vec{d} \cdot ( \vec{c} \times \vec{b} ) }{T} \\ v_0 &= \frac{ \vec{d} \cdot ( \vec{a} \times \vec{c} ) }{T} \\ w_0 &= \frac{ \vec{d} \cdot ( \vec{b} \times \vec{a} ) }{T} \\ u_x &= \frac{ c_y ( d_z - b_z ) + b_y ( c_z - d_z ) + d_y ( b_z - c_z ) }{T} \\ v_x &= \frac{ c_y ( a_z - d_z ) + a_y ( d_z - c_z ) + d_y ( c_z - a_z ) }{T} \\ w_x &= \frac{ a_y ( b_z - d_z ) + b_y ( d_z - a_z ) + d_y ( a_z - b_z ) }{T} \\ u_y &= \frac{ c_x ( b_z - d_z ) + b_x ( d_z - c_z ) + d_x ( c_z - b_z ) }{T} \\ v_y &= \frac{ c_x ( d_z - a_z ) + a_x ( c_z - d_z ) + d_x ( a_z - c_z ) }{T} \\ w_y &= \frac{ a_x ( d_z - b_z ) + b_x ( a_z - d_z ) + d_x ( b_z - a_z ) }{T} \\ u_z &= \frac{ c_x ( d_y - b_y ) + b_x ( c_y - d_y ) + d_x ( b_y - c_y ) }{T} \\ v_z &= \frac{ c_x ( a_y - d_y ) + a_x ( d_y - c_y ) + d_x ( c_y - a_y ) }{T} \\ w_z &= \frac{ a_x ( b_y - d_y ) + b_x ( d_y - a_y ) + d_x ( a_y - b_y ) }{T} \\ \end{aligned}$$ There are probably nice vector expressions for all twelve coefficients, I just have never bothered to find them out.
If this is for numerical calculation in a computer program, I would probably go with the direct Cartesian coordinate integral over the axis-aligned bounding box containing the tetrahedron, simply excluding points that are outside the tetrahedron. Note that when the face plane normals face outwards, point $\vec{p}$ is within the tetrahedron if an only if all following are true: $$\begin{aligned} \vec{n}_{abc} \cdot \vec{p} &\le d_{abc} \\ \vec{n}_{abd} \cdot \vec{p} &\le d_{abd} \\ \vec{n}_{acd} \cdot \vec{p} &\le d_{acd} \\ \vec{n}_{bcd} \cdot \vec{p} &\le d_{bcd} \\ \end{aligned}$$ If any of these are untrue, point $\vec{p}$ is not within the tetrahedron. The fraction of points inside the axis-aligned bounding box is relatively small, less than $1/6 \approx 0.16667$, but the exclusion test is very fast. The idea is that the exclusion test is so fast, that doing a lot of them is faster than doing coordinate conversions for just the points inside the tetrahedron.
In a peer-reviewed article, the volume integral over each simplex/tetrahedron is not important as it is well known, but if proof is required, the barycentric coordinate transform is straightforward to derive.