Accurate Volume Calculation Driven by Delaunay Triangulation for Coal Measurement

Volume calculation from 3D point cloud is widely used in engineering and applications. ,e existing methods either have large errors or are time-consuming. ,is paper focuses on the coal measurement. Based on the triangular mesh generated from the point cloud, each triangle is projected downward to the base plane to form a voxel. We derive the calculation formula of voxel by an integral method, which is more efficient than the method of decomposing voxel into tetrahedrons and more accurate than slicing methods. Furthermore, this paper proposes a Delaunay triangulation-driven volume calculation (DTVC) method. DTVC does not preserve the Delaunay triangles but directly calculates the volume in the process of triangulation. It saves memory and running time. Experimental results show that DTVC has achieved a good balance between error and efficiency.


Introduction
Volume calculation is to compute the space occupied by an object through point cloud that generated from 3D reconstruction. Volume calculation underpins many crucial applications, such as the volume measurement of coal [1], ore [2], earthwork [3], and tree [4,5] and the disease diagnosis in medical field [6,7]. is paper aims to the coal measurement, where volume calculation has been a challenging task due to the difficulty in determining the object boundary from the scatters of point cloud.
A rough method to calculate object volume is to fit a standard geometry, such as ellipsoid [8] and cylinder [9][10][11].
ese methods are fast and simple because they only need to fit the geometry parameters. However, due to the difference between the actual object and the regular geometry, the calculated volume error is large. e more accurate calculation methods can be roughly divided into three categories. (1) Convex hull method [12,13]. Irregular object is represented by convex hull. e convex hull can be decomposed into two triangular meshes. e corresponding projection volume is calculated by the orthographic projection method, and the difference is the object volume. e method is suitable for the convex model, and the error of nonconvex model is large. (2) Slicing method [14,15]. is method slices the point cloud using a plane that is perpendicular to a coordinate axis. e total volume was obtained by accumulating the volume of each slice. e thickness of the slice has a great influence on the result. (3) Projection method [16]. is method has two stages. First, the 2D projection of 3D point cloud is meshed, usually using the triangulation grid. Each triangle in the mesh projects to a plane. e projected triangle and the original triangle form a pentahedron. Second, the pentahedron volume is calculated, and the total volume is obtained by accumulating the volume of all the pentahedrons.
Coal pile is not a convex model, as shown in Figure 1, so the convex hull method cannot be used in this application. e top slice of coal pile has many independent areas. It is difficult to split them, so the error is large when the slicing method is used. us, the triangulation projection method is suitable for coal measurement.
ere are two methods in the pentahedron volume calculation in triangulation projection method. (1) Connect the diagonal lines of each pentahedron face, and then, the pentahedron is divided into three tetrahedrons. e volume of pentahedron is obtained by accumulating the volumes of tetrahedrons [17]. (2) e pentahedron is approximated to a triangular prism to calculate the volume [18]. ese methods have large error or are time-consuming. is paper proposes a Delaunay triangulation-driven volume calculation (DTVC) method, where the two stages in the triangulation projection method are merged together. With the construction of triangular mesh, the volume calculation is completed. e pentahedron volume calculation is also improved. It is more accurate and efficient than other methods.
We summarize our contributions as follows: (i) An accurate and direct method is proposed to calculate the volume of the pentahedron, which provides the analytic solutions of pentahedron using the integral method. (ii) We also combined the two stages to one for the triangulation projection method. e volume calculation is driven by the triangulation. (iii) e calculation error of the whole point cloud volume is less than the slicing method, and the efficiency is higher than the existing projection method. e proposed method has achieved a good balance between error and efficiency. e remainder of this paper is organized as follows. Section 2 discusses the related work about the volume calculation of point cloud. In Section 3, the accurate and direct pentahedron volume calculation is presented. In Section 4, the triangulation-driven volume calculation is introduced. Section 5 provides the experimental results. Section 6 concludes this paper.

Related Work
is section discusses the volume calculation methods.

Rough Methods.
Guo et al. [4] calculate the standing tree volume using an empirical formula, where the volume is the product of the cross-sectional area at breast height, the height plus three, and a regulatory factor. Wu et al. [19] and Okinda et al. [8] regard the jujubes and eggs as ellipsoids to calculate the volume. In [11], each mass of food on the dish is looked as a cylinder. In [9,10], the tank is approximated to a cylinder. Hadhazi et al. [7] place 16 simulated small elliptical nodules into an anthropomorphic chest phantom to calculate the lung nodule volume. Because the real objects are not standard geometry, the error is very large for these methods.

Convex Hull Method.
In 1972, Graham proposed the convex hull method [12]. Speakman and Averkov [20] introduce the convex hull volume formula of a trilinear monomial. Wu et al. [21] discuss the properties and computation of the boundary of convex hulls of cospherical circles. Alshamrani et al. [22] filter all points inside a fourvertex polygon to reduce the computational cost for computing a convex hull for a large set of points. Zhao et al. [23] introduce an incremental convex hull calculation and a fast replacement for nondominated sorting. Auat Cheein and Guivant [24] split the 3D raw data into subsets, which have the same convex hull. Calculating the convex hull from subsets can decrease the computational cost of the entire process. Kim et al. [25] propose a precise convex hull computation method for free form models, where a Gauss map organized in a hierarchy of normal pyramids and a Coons bounding volume hierarchy are prebuilt.
In application, to calculate the object volume on a conveyor belt, Song and Ou [26] apply a convex hull algorithm to find the set of the smallest convex polygon that contains all the object point cloud. Furthermore, they find a rectangle which has the smallest area to determine the object size. Bandi et al. [27] also use the convex hull method to calculate the volume of dense point cloud. Xu and Xu [13] use the incremental algorithm to calculate the convex hull of the point cloud to approximate the object. Shah et al. [28] apply convex hull to calculate the representative element of volume by simultaneously considering the two main macroscopic petrophysical parameters: porosity and singlephase permeability.

Slicing Method.
Chang et al. [29] cut the point cloud into slices of equal thickness along the z-axis. en, divide each slice along the y-axis equally and cut each slice into subintervals along the x-axis. Using the minimum and maximum coordinates in each subinterval to calculate the area of each slice. en, the slice areas are integrated along the z-axis to estimate the volume of the object. Li et al. [30] take a frame in the structured light vision as a slice. Heckel et al. [31] propose a segmentation-based partial volume analysis method to measure the volume of solid. Avdal et al. [32] estimate the flow volume by the maximum velocity in each voxel. Lyra et al. [33] add transaxial slices to measure the thyroid volume of hyperthyroid patients. Yu et al. [34] calculate the difference in the current depth image and reference depth image by splitting planes to estimate the volume of chest wall.
In [35], a slice is called a voxel, which represents a particular space that has been divided into a grid of identically sized and equally spaced cubes. Basher et al. [6] split the hippocampus into voxels. A convolutional neural network (CNN) model is designed to predict the number of voxels attributed to the hippocampus, and the number of estimated hippocampal voxels is multiplied by the voxel volume to measure the discrete volume of the hippocampus. Gao et al. [36] use the voxel-based partial volume analysis method to calculate the volume. To simulate the voxel volumes of fibres, Aronsson [37] defines each fibre as a spline curve with an elliptical cross-sectional shape and a constant twist per length unit.

Projection Method.
Potena et al. [38] adopt a tailored meshing strategy that performs a leaf-based clustering to obtain an approximated canopy estimation and an iterative cluster connection. e volume is estimated by the resulting mesh. In volume estimation of nongeometric shape cavity, Toumpaniaris et al. [39] design the distances for the triangulation are from a random point of the cavity to the inner wall.
e rand point and the triangle form a triangular pyramid. e volume summation of the triangular pyramids is the volume of the irregular body. Cai et al. [40] measure the potato volume using laser triangulation in cylindrical coordinates.
Hu et al. [17] construct the mesh of 3D point cloud by using implicit B-spline. e triangles of each facet are projected to the xOy plane, and the projected triangles are combined with the corresponding facet triangles to form a pentahedron. e pentahedron is divided into three tetrahedrons by connecting the diagonals of each surface. en, the pentahedron volume is computed by summarizing the tetrahedron volumes. In [16], the 3D liver model which is constructed by the marching cube algorithm is composed of a number of triangle meshes, and the normal vectors of different directions of all triangle meshes are modified, traversed, and projected onto a projection plane of 3D spatial coordinates one by one. en, a lot of pentahedrons can be built, and the volumes of all pentahedrons are calculated by the subdivision method, and the final volume of 3D liver model is the algebraic sum of all pentahedron volumes. Wang et al. [18] transform the pentahedron approximately into a triangular prism to calculate the volume. Taking the triangle as the smallest unit, the space irregular body is projected to the xOy plane, and the area of the projected triangle is calculated by using the Helen formula. e volume of the Mitsubishi column can be obtained by taking the average value of the Z coordinate value of the projection point as the height.

Accurate Volume Calculation Method
is section describes the detail of accurate volume calculation method, which is based on the triangulation projection method.

Volume Calculation Method.
e coal pile is stacked on the horizontal ground. e ground is called the base plane. e 3D point cloud generated from the laser scanner or camera array is the set of discrete surface points. Using the triangulation method, the messy surface points can be connected to form a triangular mesh, as shown in Figure 2.
We choose a triangle △ABC in the triangular mesh and project it to the base plane to get △A ′ B ′ C ′ . e triangles △ABC and △A ′ B ′ C ′ and the quadrilaterals ABB ′ A ′ , BCC ′ B ′ , and ACC ′ A ′ form a pentahedron, which is called voxel in this paper.
Suppose the voxel set is σ i n i�1 and the volume of σ i is v i , and then, the coal volume can be computed as follows:

Equation of the Top Surface.
e equation of the top surface of σ i is required before calculating the volume using integral method. Take the base plane as the xOy plane and zaxis pointing upward to establish the coordinate system. We rearrange the order of A i , B i , and C i so that x i1 ≤ x i2 ≤ x i3 . From the view of the positive direction of zaxis, B i can be on the right or left side of A i C i , as shown in Figure 3. We can calculate the normal vector of the top plane, as shown in equation (2), and is abbreviated as in equation (3): where which can be transformed as follows:

Scientific
Programming e coordinates of A i ′ , B i ′ , and C i ′ are (x i1 , y i1 , 0), (x i2 , y i2 , 0), and (x i3 , y i3 , 0), respectively. We first only consider the case of x i1 < x i2 < x i3 , and then, the equations of A i ′ B i ′ , B i ′ C i ′ , and A i ′ C i ′ can be computed by equations (6), (7), and (8), respectively:

e Edge
where where where

e Volume Calculation of Voxel.
After the equations of the top surface and the projected triangle edges are calculated, the volume of voxel σ i can be computed using integral method. As shown in Figure 4, the projected triangle is divided into two parts by the blue line L. It is corresponding that the pentahedron is divided into two parts by a plane that is perpendicular to the xOy plane and passes through the line L. e two parts can be calculated by equations (9) and (10), respectively. en, we can compute the volume of voxel σ i , which is shown in equation (11): Let δ � k 3 x i2 + b 3 , as the red part that is shown in Figure 4. By solving these integral equations, we can get the following equations: In Section 3.3, we limit the condition of

Triangulation-Driven Volume Calculation
is section introduces the DTVC algorithm. DTVC finds triangles one by one without repetition in the 3D point cloud. With the discovery of the triangle, the volume of its projected pentahedron is calculated at the same time. In this algorithm, we need not to create a triangular mesh and a very large storage structure of triangle set, which can save storage space and time in the actual scene. [41] is a popular used method to construct triangulation mesh. Suppose the point set is V, and T is any triangle of V. If T is a Delaunay triangle of V, only the interior of the circumcircle of each triangle in T does not contain any point in V, which is called the empty circumcircle principle. e point-by-point insertion algorithm [42,43], also called the Bowyer-Watson algorithm, is the most widely used Delaunay triangulation algorithm.

Delaunay Triangulation. Delaunay algorithm
Firstly, a large triangle is built to contain all the points, which is called the supertriangle. en, insert a point into the triangular, which is connected with the three vertices of the triangle to form three new triangles. ey are checked one by one according to the empty circumcircle principle. Repeat these steps until all points are checked.  Scientific Programming

Supertriangle Creation.
e supertriangular is a special triangle which contains all the points in V.
is section introduces the creation of the supertriangular.
Firstly, find the maximum and minimum coordinates of all points in x and y directions, which are denoted as x max , x min , y max , and y min .
ese four values construct a rectangle, which contains all the points in V. e rectangle is the minimum bounding box of all points, which is shown as green rectangle in Figure 5.
Secondly, the bounding box is inflated for ϵ length on four directions to form the red rectangle ABC D. According to equation (14), the coordinates of A, B, C, and D are (u, t), (v, t), (v, s), and (u, s), respectively. Some points may be on the edge of the green rectangle, but all points are inside of the red rectangle: where ε is a very small positive number.
Fourthly, connect and extend CA and FB, which intersect at point G.
en, the coordinate of G is ((u + v/2), 2(t − s)).  Figure 5: e supertriangle. e green rectangle is the minimum bounding box of all points. e red rectangle is obtained by the bounding box inflating for ϵ length on four directions. e blue triangle is the supertriangle. Scientific Programming 5 e blue triangle △EFG is the supertriangle.

Triangulation-Driven Algorithm.
e key of triangulation is to determine whether a triangle is a Delaunay triangle. In Figure 6(a), △AB D, △AC D, and △BC D are undetermined triangles. When inserting a point E, we make circumcircle for each undetermined triangle, as shown in Figure 6(b). According to the relative position of E and circumcircle, we can determine the triangle type. e decision rules are as follows: In the DTVC algorithm, instead of storing the Delaunay triangles, we directly calculate the projected volume. e algorithm is shown in Algorithm 1. is algorithm is used to find triangle without repetition. We do not create and store the mesh in the calculation. When all the points in the 3D point cloud are traversed, the volume is obtained.
ere is a shared algorithm to calculate volume of a set, which is shown in Algorithm 1. e calculation of this algorithm is based on equations (12) and (13).

Datasets.
In order to analyze the error of the algorithm, we simulate two point cloud data, as shown in Figure 7. To test the volume calculation in real scene, we take 59 photos in a 320 × 110 × 49 coal shed to generate the point cloud. e resolution of the photos is 3072 × 2048.

Experiment Setup.
In error analysis, we use Matlab R2018b implementation tool. e processor and memory are Intel Core i7-8565U and 16.0 GB, respectively. e operating system is Win10.
In the real scene, the 3D reconstruction and volume calculation system run on Ubuntu 18.04 LTS. e CPU, GPU, and memory are Xeon Gold 6226R, Nvidia 16 GB Tesla T4, and 256G DDR4, respectively. e error is calculated by equation (15), and the CUP time is the difference between the finish time and the start time of the algorithm. We repeat each algorithm 100 times and average the errors and CUP times as the final experiment result: where v and v r are the measured and the real volume, respectively.

Compare with State-of-the-Art Method.
We choose some popular used methods to compare with our proposed method. e notation for these methods is as follows: (i) Slice-# [29]. Slicing method, and "#" is the slice number. (ii) P-T3 [17]. e projection method which divides the pentahedron into 3 tetrahedrons. (iii) DTVC. Our Delaunay triangulation-driven volume calculation.
Tables 1 and 2 show the result of single peak and multiple peaks, respectively. Experimental results show that the slicing method has the most efficient operation. However, the error rate is also large. When the slice number is increased, the error goes down and the running time rises up. e results based on the slicing model are difficult to be applied to practical solutions. P-T3 and our DTVC have the same low error, but P-T3 is more time-consuming than DTVC because P-T3 requires three determinant calculations of matrix. P-T3 spends more than 2 times on the running time evaluation.
Experimental results show that our DTVC method has achieved a good balance between error and efficiency.

Application in Real Scene.
In a real scene, we deploy a camera matrix in a coal shed. We generate the point cloud using the 3D movement structure recovery method. Based on the point cloud, we use the proposed DTVC method to calculate the volume of coal pile. Figure 8 shows a display interface of the system.

Conflicts of Interest
e authors declare that they have no conflicts of interest.