Fast volume computing algorithm of arbitrary polyhedron in 3D space
Introduction
Computing of the volume of polyhedron in 3D space isn’t a trivial problem, but the following trivial method exists: dividing the volume into simple pyramids and counting the sum of these volumes. However, this methodology is difficult for implementation, as well as it is resource-intensive and slow. What is more, computing algorithm according to this methodology is difficult to parallelize. And taking into account trends in the field of parallel computing on GPU development, the ability to increase the rate of computing repeatedly is being lost.
In this article the algorithm for determining the volume of arbitrary geometry in 3D space in terms of fast computing is described.
Input data
Preparatory step. A rectangular grid (pos.1 on fig.1) around geometry (pos.2 on fig.1) is built using its minimum and maximum coordinates.

User chooses rectangular grid parameters. For example, the number of grid cells can be taken equal to the number of stream multi-processors on GPU. Therefore, our task turns into computing and subsequent summing of the volumes occupied by fragments of geometries in cells, as shown on Figure 2.

Practically in all cases object geometry is a corresponding set of triangles with an appropriate bypass (or normal). This set can be stored as a whole massive and you may appeal to it parallel from every cell or perform geometry pre-processing for all cells. In this example we don’t describe geometry pre-processing for cells, because this is the subject of a separate topic for optimization.
Algorithm of cells volume definition:
Step 1. Carry out mesh clipping by means of box (grid cell). To simplify and accelerate the algorithm, clipping is performed with the help of 6 planes that form a box. Clipping is performed by each plane alternately in any order.
Step 2. To project onto XOY all polygons which make clipped mesh. Compute the volume under each triangle. The volumes are computed as the sums of three pyramid volumes, which form the shape (see Fig. 3).

Step 3. Depending on the orientation of normal these partial volumes are included into the sum with opposite signifier (signifier of scalar product of the normal of polygon and projection plane).

Step 4. Compute cross-sectional area formed by upper face box (parallel to projection plane) and mesh (see figure 5).

In order to compute the entire volume of the object inside the box it is necessary not only to cut it by means of the box, but also to close it (which is a very time consuming task). In order to gain the volume all rectangles are projected onto the box face. Therefore, it is possible to simplify circuit computing, having processed only the upper face, because projections of other faces on box face give zero area.
Compute cross-sectional area (see Fig. 6) with modification described below.

Modification takes into account unclosed cases (polygon can cross edges of the upper face box, in this case it will cross faces of rectangle) (see fig. 7).

In such case the upper edge of the face should form two-dimensional polygon locking. Again, there is considered the fact that it is no need to compute locking on all edges of the upper face, it will be enough to compute locking only on the upper edge. Thus, next operations take place on one segment. Further, the term edge will mean namely edge A on Fig. 8.

As current subproblem has been reduced to one-dimensional, then all geometric concepts are used in one-dimensional sense (projection on the line). For example, the normal (or any other normalized vector), in fact, is 1 or -1. In connection with such binary structure normals can be divided into left (1) and right (+1). Figure 9 illustrates the scheme of computing the locking on top edge.

The edge is divided into segments and each segment is determined by the number of normals. We mean that the number of normals of one segment is the number of normals directed from the given segment minus the number of normals directed to this segment, it is also necessary to take into account only those normals of mesh triangles, which have exactly one intercut with the edge.
Initially, the whole edge forms a single source segment with the number of normals equal to 0. Depending on the parity number of triangles, which cross the edge, internal (relative to the mesh) edge segments may have number of normals equal to 0 or 1. Depending on the implementation of the previous steps of the algorithm some boundary situations may occur, for example, as on Figure 10:

Let us consider the case when triangle intersects with the edge on a segment AB (more than one point of intersection). Segment AB needs to be locked from further computing of normals for two reasons:
• The area under this segment is taken into account because this segment is on the list of two-dimensional faces of the polygon;
• It is impossible to project the normal onto the edge.
Algorithm completion.
On the first step of the algorithm on determination of mesh cross-section with each cell, there are cases when cross-sections of the mesh with the cell do not exist (the cell is either inside or outside). If a cross-section of mesh with a cell exists, these cells are marked as empty (free).
Further, it is necessary to mark empty cells as the cells existing either in the volume or out of the volume. For example, to mark with 0 if the cells are out of the volume or with the value of cell volume, if the amount is included in the value.
You can use one of the two algorithms or their combination:
• Row-wise scan algorithm;
• Algorithm of recursive volume fill.
Both of these algorithms are widely discussed on the Internet. The essence of the first algorithm is the following: arbitrary direction is selected (eg, along OZ axis).The plane along the chosen direction is consistently chosen (for example, XOY plane in case of OZ axis). You begin to scan on the plane each line which is a consistent set of cells along a particular direction. Cells marking is reduced to highlighting empty cells until any filled cell is met. This algorithm has many restrictions that need to be considered separately (if only 1 cell on the line is filled; if 2 cells are filled, but they are top of geometry, etc.).
The essence of the second algorithm is reduced to the choice of any empty cell and then to the recursive analysis of neighbor cells and their subsequent addition to the sampling queue. If neighbor cell is full, this cell is not included in the queue for the sample. As soon as it is impossible to choose neighbor cell, the whole sequence of cells is marked as selected cells in a certain volume. Next arbitrary cell is selected and all the same steps are conducted. As a result, empty cells are marked as belonging to a particular volume. Further it is necessary to determine whether the selected volume is empty or filled. This is the complexity of the algorithm. The simplest case is the following one: in the course of recursion as soon as the first full cell or volume border is there, the analysis of the entire volume inside or outside is performed. Analysis is performed on the basis of the normals of triangles that make up the mesh in an occupied cell. But this operation is either time-consuming or resource-intensive. Due to the drawbacks of both methodologies, it is advisable to resort to their combination. For example, to use only a recursive volume marking and then pass the cell line with marked volumes. Such technology can also be quite easily parallelized.
Conclusion.
Thus, the algorithm of computing the volume of any multi-segment geometry in space has been introduced. This algorithm can be used for determination of oil or gas amount in the field. Apparent advantage of the algorithm is its ability to parallelization and as a result – rapid result achievement. The disadvantage of the algorithm is complexity of its understanding and implementation.
It would be great to discuss practical aspects of the implementation in different architetures, like GPUs for example...
Vicente's comment from LinkedIn: "very nice article, great explanation...it woud be interesting to discuss practical impementation in gpus for example.
Although new nvidia devices allow sort of recursive algorithms, most deveices require a more explicit approach."
Hello Vicente!
We intend to apply this algorithm on multi-core processors either on GPUs. Moreover for GPGPU application we usually use Nvidia adapters. Currently we actively test this algorithm; GPU architecture consideration will be later.
Computing the volume of a polyhedron is a trivial task assuming you have a solid: create a tessellation of the boundary, calculate the volume of the tetrahedron defined by each triangle on the boundary and a fix point in space, add all the volumes, and you are done. This method is very efficient even when you have millions of triangles. I do that all the time. I don't see a need for anything more complicated. If you are concerned about errors there are techniques to improve the accuracy of the summation of the volumes of each tetrahedron.
The algorithm that you have already described is post factum "classic" algorithm of calculating the volume of a polyhedron, which is true only for a single object (polyhedron).
Comparison between the proposed "Fast volume computing algorithm "and "classic volume computing algorithm" is the following:
The proposed Fast volume computing algorithm can compute the volume, occupied by a number of intersected polyhedra with the inaccuracy only in those cells in which there is contact of polyhedra faces.
The Opportunities for parallelization of the algorithm is not worse the opportunities of the classic algorithm that you have already described.
Due to the fact that the algorithm was written for solving a particular objective (in our objective it is necessary to operate with the resulting mesh and the volumes occupied by the objects in the cell), then many researchers may be interested in the following: 1) the method of calculating the volume of the cell, not inferior to the "classic" method either in speed or in accuracy, and 2) a considerable amount of "related" information obtained in the course of the algorithm (eg, mesh, triangles cut off by grid cell, etc.).
With all due respect, this article loosely describes a method to approximate the volume of a polyhedron using voxels; there is no mention on how the technique works on intersecting polyhedra. (If the polyhedra are mutually disjoint what you call classic volume calculation is still valid, and will be much faster than the method proposed here.)
There are no error estimates associated to the error induced by the volume calculation of the cells near the boundary. If the error estimate is simply the sum of the volumes of the cells that intersect the boundary, we are considering a relatively large error.
I cannot think of any application where the volume of intersecting polyhedra need to be calculated without actually calculating the actual intersection or union of the polyhedra. Is there a specific application for the method being proposed?
The method described in this article, is used to compute the volume of body intersection in dynamics. You can build three-dimensional cell representations of the simulated objects and then intersect them (deform). Of course, you will get the error, but the error depends on the cells size. This algorithm provides accuracy-speed tradoff with the help of regulation of only one parameter - the cells size. The presence of the "real" amount of cell body allows to use such representation in modeling systems.