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.

  • 1
    $\begingroup$ @Adag1242: Thanks for the edit; I had bugs there, but you fixed them. I noticed additional bugs, missing $u$ factors in the barycentric coordinate calculation. $\endgroup$ – Anonymous Coward Apr 20 '20 at 15:34
  • 1
    $\begingroup$ @Adag1242: No, that is not a stupic question. The problem is that the $(0,0,0)$,$(1,0,0)$,$(0,1,0)$,$(0,0,1)$ tetrahedron is likely skewed with respect to the actual tetrahedron. This is what the $V_d$ factor (which is the determinant of the partial derivatives matrix) "corrects", and it does depend on the coordinates of the tetrahedral vertices. So, I don't think what you suggest can be used to find the exact solution. $\endgroup$ – Anonymous Coward Apr 20 '20 at 15:41
  • 1
    $\begingroup$ Do note that if you have a function $$\delta(x, y, z) = \begin{cases} 1, & \text{inside the tetrahedron} \\ 0, & \text{elsewhere} \end{cases}$$you can write the volume integral as$$\int_{z_{min}}^{z_{max}} \int_{y_{min}}^{y_{max}} \int_{x_{min}}^{x_{max}} \delta(x, y, z) g(x, y, z) ~ d x ~ d y ~ d z$$where the limits specify an axis-aligned cuboid that completely contains the tetrahedron. If the cuboid is larger, you only waste a bit of computational power. If you store the tetrahedral face normal vectors and signed distances, $\delta(x, y, z)$ is easy to implement. $\endgroup$ – Anonymous Coward Apr 20 '20 at 15:48
  • 1
    $\begingroup$ @Adag1242: Yes, except barycentric, not Bayesian (unless there is something important you're not telling us). Note that $$V_d(u, v, w) = V_{vww} v w^2 + V_{vw} v w + V_{w} w + V_{v} v + V_0$$which is part of the integral, and I didn't show in explicit form; the constants $V_{*}$ depend on the vertex coordinates only. If you evaluate numerical integrals only, then using the exclusion method ($\delta(x,y,z)$) and Monte Carlo integration might be useful. $\endgroup$ – Anonymous Coward Apr 20 '20 at 16:09
  • 1
    $\begingroup$ (At the risk of derailing your idea, I would also consider spatial subdivision, instead of tetrahedral mesh. The idea is that you recursively split the cuboid volume in half, along its longest axis (along one of the coordinate axes), until the average within the cuboid sufficiently describes the function within it. Note that for numerical solutions, you can uniform-randomly sample $g(x,y,z)$ within the cuboid to determine whether subdivision is necessary (based on variance, or just min/max), before calculating the integral at all.) $\endgroup$ – Anonymous Coward Apr 20 '20 at 16:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.