# Accurate Volume Calculation Driven by Delaunay Triangulation for Coal Measurement

Authors: Yong Liu; Yanwei Zheng
Journal: Scientific Programming
 (2021)

---

## Abstract

Volume calculation from 3D point cloud is widely used in engineering and applications. The existing methods either have large errors or are time-consuming. This 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.


---

## Body

## 1. 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]. This 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–11]. These 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.The more accurate calculation methods can be roughly divided into three categories. (1) Convex hull method [12, 13]. Irregular object is represented by convex hull. The convex hull can be decomposed into two triangular meshes. The corresponding projection volume is calculated by the orthographic projection method, and the difference is the object volume. The method is suitable for the convex model, and the error of nonconvex model is large. (2) Slicing method [14, 15]. This method slices the point cloud using a plane that is perpendicular to a coordinate axis. The total volume was obtained by accumulating the volume of each slice. The thickness of the slice has a great influence on the result. (3) Projection method [16]. This 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. The 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 Figure1, so the convex hull method cannot be used in this application. The 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. Thus, the triangulation projection method is suitable for coal measurement.Figure 1
Coal pile from different view positions. It is not a convex model.There 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. The volume of pentahedron is obtained by accumulating the volumes of tetrahedrons [17]. (2) The pentahedron is approximated to a triangular prism to calculate the volume [18]. These methods have large error or are time-consuming.This 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. The 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. The volume calculation is driven by the triangulation.(iii)
The calculation error of the whole point cloud volume is less than the slicing method, and the efficiency is higher than the existing projection method. The proposed method has achieved a good balance between error and efficiency.The remainder of this paper is organized as follows. Section2 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.


## 2. Related Work

This section discusses the volume calculation methods.


### 2.1. 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.


### 2.2. 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 four-vertex 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 single-phase permeability.


### 2.3. Slicing Method

Chang et al. [29] cut the point cloud into slices of equal thickness along the z-axis. Then, 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. Then, 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.


### 2.4. 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. The 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. The rand point and the triangle form a triangular pyramid. The 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. The 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. The pentahedron is divided into three tetrahedrons by connecting the diagonals of each surface. Then, 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. Then, 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. The 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.


## 2.1. 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.


## 2.2. 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 four-vertex 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 single-phase permeability.


## 2.3. Slicing Method

Chang et al. [29] cut the point cloud into slices of equal thickness along the z-axis. Then, 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. Then, 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.


## 2.4. 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. The 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. The rand point and the triangle form a triangular pyramid. The 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. The 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. The pentahedron is divided into three tetrahedrons by connecting the diagonals of each surface. Then, 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. Then, 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. The 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.


## 3. Accurate Volume Calculation Method

This section describes the detail of accurate volume calculation method, which is based on the triangulation projection method.


### 3.1. Volume Calculation Method

The coal pile is stacked on the horizontal ground. The ground is called the base plane. The 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 Figure2.Figure 2
The triangular mesh and voxel. Triangular mesh is generated after triangulation of point cloud. The triangle△A′B′C′ on the base plane is the projection of original triangle △ABC in the triangular mesh. The pentahedron ABC−A′B′C′ is called a voxel.We choose a triangle△ABC in the triangular mesh and project it to the base plane to get △A′B′C′. The 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σii=1n and the volume of σi is vi, and then, the coal volume can be computed as follows:(1)V=∑i=1nvi.


### 3.2. Equation of the Top Surface

The 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 z-axis pointing upward to establish the coordinate system. Suppose the top and bottom surfaces of voxel σi are △AiBiCi and △Ai′Bi′Ci′, respectively. The coordinates of Ai, Bi, and Ci are xi1,yi1,zi1, xi2,yi2,zi2, and xi3,yi3,zi3, respectively.We rearrange the order ofAi, Bi, and Ci so that xi1≤xi2≤xi3. From the view of the positive direction of z-axis, Bi can be on the right or left side of AiCi, 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):(2)ni⟶=i⟶j⟶k⟶xi2−xi1yi2−yi1zi2−zi1xi3−xi1yi3−yi1zi3−zi1,(3)ni⟶=a,b,c,where a=yi2−yi1zi3−zi1−yi3−yi1zi2−zi1, b=xi3−xi1zi2−zi1−xi2−xi1zi3−zi1, and c=xi2−xi1yi3−yi1−xi3−xi1yi2−yi1.Figure 3
Rearrange the order ofAi, Bi, and Ci so that xi1≤xi2≤xi3, where Ai, Bi, and Ci are the vertices of voxel σi and xi1,yi1,zi1,xi2,yi2,zi2,xi3,yi3,zi3 are the coordinates of Ai,Bi,Ci, respectively. (a) Bi′ is on the right of the line Ai′Ci′ in the base plane. (b) Bi′ is on the left of the line Ai′Ci′ in the base plane.
(a)(b)In△Ai′Bi′Ci′, Ai′, Bi′, and Ci′ are not collinear, so Ai, Bi, and Ci are not in a plane that is perpendicular to the base plane. Thus, c≠0.The equation of the top plane can be expressed as(4)ax−xi1+by−yi1+cz−zi1=0,which can be transformed as follows:(5)z=a′x+b′y+c′,where a′=−a/c, b′=−b/c, and c′=a/cxi1+b/cyi1+zi1.


### 3.3. The Edge Equation of the Projected Triangle

Looking from the positive to negative direction ofz-axis, Bi′ can be on the right or left of line Ai′Ci′, as shown in Figures 3(a) and 3(b), respectively. They are corresponding to the projections, as shown in Figures 4(a) and 4(b), respectively.Figure 4
The projection in the base plane. (a)Bi′ is above the line Ai′Ci′. (b) Bi′ is under the line Ai′Ci′.
(a)(b)The coordinates ofAi′, Bi′, and Ci′ are xi1,yi1,0, xi2,yi2,0, and xi3,yi3,0, respectively. We first only consider the case of xi1<xi2<xi3, and then, the equations of Ai′Bi′, Bi′Ci′, and Ai′Ci′ can be computed by equations (6), (7), and (8), respectively:(6)y=k1x+b1,where k1=yi2−yi1/xi2−xi1 and b1=xi2yi1−xi1yi2/xi2−xi1.(7)y=k2x+b2,where k2=yi3−yi2/xi3−xi2 and b2=xi3yi2−xi2yi3/xi3−xi2.(8)y=k3x+b3,where k3=yi3−yi1/xi3−xi1 and b3=xi3yi1−xi1yi3/xi3−xi1.


### 3.4. The 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 Figure4, 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. The two parts can be calculated by equations (9) and (10), respectively. Then, we can compute the volume of voxel σi, which is shown in equation (11):(9)vi1=∫xi1xi2dx∫k3x+b3k1x+b1dy∫0a′x+b′y+c′dz,(10)vi2=∫xi2xi3dx∫k3x+b3k2x+b2dy∫0a′x+b′y+c′dz,(11)vi=vi1+vi2.Letδ=k3xi2+b3, as the red part that is shown in Figure 4. By solving these integral equations, we can get the following equations:(12)vi1=a′yi2−δ36k2−k32−a′xi2+c′yi2−δ22k2−k3+b′yi33−yi236k2−b′yi33−δ36k3,(13)vi2=−a′yi2−δ36k1−k32+a′xi2+c′yi2−δ22k1−k3+b′yi33−yi136k1+b′yi13−δ36k3.In Section3.3, we limit the condition of xi1<xi2<xi3. If xi1=xi2 or xi2=xi3, then there is vi1=0 or vi2=0, respectively.


## 3.1. Volume Calculation Method

The coal pile is stacked on the horizontal ground. The ground is called the base plane. The 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 Figure2.Figure 2
The triangular mesh and voxel. Triangular mesh is generated after triangulation of point cloud. The triangle△A′B′C′ on the base plane is the projection of original triangle △ABC in the triangular mesh. The pentahedron ABC−A′B′C′ is called a voxel.We choose a triangle△ABC in the triangular mesh and project it to the base plane to get △A′B′C′. The 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σii=1n and the volume of σi is vi, and then, the coal volume can be computed as follows:(1)V=∑i=1nvi.


## 3.2. Equation of the Top Surface

The 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 z-axis pointing upward to establish the coordinate system. Suppose the top and bottom surfaces of voxel σi are △AiBiCi and △Ai′Bi′Ci′, respectively. The coordinates of Ai, Bi, and Ci are xi1,yi1,zi1, xi2,yi2,zi2, and xi3,yi3,zi3, respectively.We rearrange the order ofAi, Bi, and Ci so that xi1≤xi2≤xi3. From the view of the positive direction of z-axis, Bi can be on the right or left side of AiCi, 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):(2)ni⟶=i⟶j⟶k⟶xi2−xi1yi2−yi1zi2−zi1xi3−xi1yi3−yi1zi3−zi1,(3)ni⟶=a,b,c,where a=yi2−yi1zi3−zi1−yi3−yi1zi2−zi1, b=xi3−xi1zi2−zi1−xi2−xi1zi3−zi1, and c=xi2−xi1yi3−yi1−xi3−xi1yi2−yi1.Figure 3
Rearrange the order ofAi, Bi, and Ci so that xi1≤xi2≤xi3, where Ai, Bi, and Ci are the vertices of voxel σi and xi1,yi1,zi1,xi2,yi2,zi2,xi3,yi3,zi3 are the coordinates of Ai,Bi,Ci, respectively. (a) Bi′ is on the right of the line Ai′Ci′ in the base plane. (b) Bi′ is on the left of the line Ai′Ci′ in the base plane.
(a)(b)In△Ai′Bi′Ci′, Ai′, Bi′, and Ci′ are not collinear, so Ai, Bi, and Ci are not in a plane that is perpendicular to the base plane. Thus, c≠0.The equation of the top plane can be expressed as(4)ax−xi1+by−yi1+cz−zi1=0,which can be transformed as follows:(5)z=a′x+b′y+c′,where a′=−a/c, b′=−b/c, and c′=a/cxi1+b/cyi1+zi1.


## 3.3. The Edge Equation of the Projected Triangle

Looking from the positive to negative direction ofz-axis, Bi′ can be on the right or left of line Ai′Ci′, as shown in Figures 3(a) and 3(b), respectively. They are corresponding to the projections, as shown in Figures 4(a) and 4(b), respectively.Figure 4
The projection in the base plane. (a)Bi′ is above the line Ai′Ci′. (b) Bi′ is under the line Ai′Ci′.
(a)(b)The coordinates ofAi′, Bi′, and Ci′ are xi1,yi1,0, xi2,yi2,0, and xi3,yi3,0, respectively. We first only consider the case of xi1<xi2<xi3, and then, the equations of Ai′Bi′, Bi′Ci′, and Ai′Ci′ can be computed by equations (6), (7), and (8), respectively:(6)y=k1x+b1,where k1=yi2−yi1/xi2−xi1 and b1=xi2yi1−xi1yi2/xi2−xi1.(7)y=k2x+b2,where k2=yi3−yi2/xi3−xi2 and b2=xi3yi2−xi2yi3/xi3−xi2.(8)y=k3x+b3,where k3=yi3−yi1/xi3−xi1 and b3=xi3yi1−xi1yi3/xi3−xi1.


## 3.4. The 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 Figure4, 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. The two parts can be calculated by equations (9) and (10), respectively. Then, we can compute the volume of voxel σi, which is shown in equation (11):(9)vi1=∫xi1xi2dx∫k3x+b3k1x+b1dy∫0a′x+b′y+c′dz,(10)vi2=∫xi2xi3dx∫k3x+b3k2x+b2dy∫0a′x+b′y+c′dz,(11)vi=vi1+vi2.Letδ=k3xi2+b3, as the red part that is shown in Figure 4. By solving these integral equations, we can get the following equations:(12)vi1=a′yi2−δ36k2−k32−a′xi2+c′yi2−δ22k2−k3+b′yi33−yi236k2−b′yi33−δ36k3,(13)vi2=−a′yi2−δ36k1−k32+a′xi2+c′yi2−δ22k1−k3+b′yi33−yi136k1+b′yi13−δ36k3.In Section3.3, we limit the condition of xi1<xi2<xi3. If xi1=xi2 or xi2=xi3, then there is vi1=0 or vi2=0, respectively.


## 4. Triangulation-Driven Volume Calculation

This 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.


### 4.1. Delaunay Triangulation

Delaunay algorithm [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. The point-by-point insertion algorithm [42, 43], also called the Bowyer–Watson algorithm, is the most widely used Delaunay triangulation algorithm.Firstly, a large triangle is built to contain all the points, which is called the supertriangle. Then, insert a point into the triangular, which is connected with the three vertices of the triangle to form three new triangles. They are checked one by one according to the empty circumcircle principle. Repeat these steps until all points are checked.


### 4.2. Supertriangle Creation

The supertriangular is a special triangle which contains all the points inV. This section introduces the creation of the supertriangular.Firstly, find the maximum and minimum coordinates of all points inx and y directions, which are denoted as xmax,xmin,ymax, and ymin. These four values construct a rectangle, which contains all the points in V. The rectangle is the minimum bounding box of all points, which is shown as green rectangle in Figure 5.Figure 5
The supertriangle. The green rectangle is the minimum bounding box of all points. The red rectangle is obtained by the bounding box inflating forϵ length on four directions. The blue triangle is the supertriangle.Secondly, the bounding box is inflated forϵ length on four directions to form the red rectangle ABCD. 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:(14)u=xmin−ɛ,v=xmax+ɛ,s=ymin−ɛ,t=ymax+ɛ,where ɛ is a very small positive number.Thirdly, extendCD to E, which satisfies DE=1/2CD. Extend DC to F, which satisfies CF=1/2CD. Then, the coordinates of E and F are 3u−v/2,s and 3v−u/2,s, respectively.Fourthly, connect and extendCA and FB, which intersect at point G. Then, the coordinate of G is u+v/2,2t−s.The blue triangle△EFG is the supertriangle.


### 4.3. Triangulation-Driven Algorithm

The key of triangulation is to determine whether a triangle is a Delaunay triangle. In Figure6(a), △ABD, △ACD, and △BCD 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. The decision rules are as follows:(i)
PointE is on the right of circumcircle of △BCD (green circle), and then, △BCD is a Delaunay triangle.(ii)
PointE is in the circumcircle of △ACD (blue circle), and then, △ACD is not a Delaunay triangle. Delete △ACD from the undetermined triangles, and △ADE, △ACE, and △CDE are new undetermined triangles.(iii)
PointE is outside and not on the right of the circumcircle of △ABD (red circle), and then, △ABD is still an undetermined triangle.Figure 6
Delaunay triangulation. When pointE is inserted, the undetermined triangles are checked by circumcircle testing. (a) Undetermined triangles. (b) Circumcircle test.
(a)(b)In the DTVC algorithm, instead of storing the Delaunay triangles, we directly calculate the projected volume. The algorithm is shown in Algorithm1. This 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.Algorithm 1: Delaunay triangulation-driven volume calculation.
(i)
Input: the vertex list V(ii)
Output: the volume Vol of V(1)
Vol=0;(2)
Initialize the temporary triangle listT;(3)
Calculate the supertriangle△EFG from V;(4)
V⋅AddE,V⋅AddF,V⋅AddG;(5)
T⋅Add△EFG;(6)
foreachPinVdo(7)
Initialize the edge bufferEdge;(8)
foreach△ABCinTdo(9)
Calculate the triangle circumcircle centerO and radius r;(10)
ifP⋅x>O⋅x+rthen(11)
Vol+= Algorithm 2T,E,F,G;(12)
end(13)
ifOP≤rthen(14)
Edge⋅AddAB,Edge⋅AddBC,Edge⋅AddAC;(15)
T⋅Delete△ABC;(16)
end(17)
end(18)
Delete all doubly specified edges fromEdge;(19)
foreachABinEdgedo(20)
T⋅Add△ABP;(21)
end(22)
end(23)
Vol+= Algorithm 2T,E,F,G;There is a shared algorithm to calculate volume of a set, which is shown in Algorithm1. The calculation of this algorithm is based on equations (12) and (13).Algorithm 2: Triangle projection volume calculation.
(i)
Input: the triangle list T, and the excluded point set P(ii)
Output: the summary volume Vol of T, where the triangle is excluded in which one of vertices is in P(1)
Vol=0;(2)
foreach△ABCinTdo(3)
ifA∉PandB∉PandC∉Pthen(4)
Calculate the projection volumevi of △ABC;(5)
Vol+=vi;(6)
end(7)
end


## 4.1. Delaunay Triangulation

Delaunay algorithm [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. The point-by-point insertion algorithm [42, 43], also called the Bowyer–Watson algorithm, is the most widely used Delaunay triangulation algorithm.Firstly, a large triangle is built to contain all the points, which is called the supertriangle. Then, insert a point into the triangular, which is connected with the three vertices of the triangle to form three new triangles. They are checked one by one according to the empty circumcircle principle. Repeat these steps until all points are checked.


## 4.2. Supertriangle Creation

The supertriangular is a special triangle which contains all the points inV. This section introduces the creation of the supertriangular.Firstly, find the maximum and minimum coordinates of all points inx and y directions, which are denoted as xmax,xmin,ymax, and ymin. These four values construct a rectangle, which contains all the points in V. The rectangle is the minimum bounding box of all points, which is shown as green rectangle in Figure 5.Figure 5
The supertriangle. The green rectangle is the minimum bounding box of all points. The red rectangle is obtained by the bounding box inflating forϵ length on four directions. The blue triangle is the supertriangle.Secondly, the bounding box is inflated forϵ length on four directions to form the red rectangle ABCD. 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:(14)u=xmin−ɛ,v=xmax+ɛ,s=ymin−ɛ,t=ymax+ɛ,where ɛ is a very small positive number.Thirdly, extendCD to E, which satisfies DE=1/2CD. Extend DC to F, which satisfies CF=1/2CD. Then, the coordinates of E and F are 3u−v/2,s and 3v−u/2,s, respectively.Fourthly, connect and extendCA and FB, which intersect at point G. Then, the coordinate of G is u+v/2,2t−s.The blue triangle△EFG is the supertriangle.


## 4.3. Triangulation-Driven Algorithm

The key of triangulation is to determine whether a triangle is a Delaunay triangle. In Figure6(a), △ABD, △ACD, and △BCD 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. The decision rules are as follows:(i)
PointE is on the right of circumcircle of △BCD (green circle), and then, △BCD is a Delaunay triangle.(ii)
PointE is in the circumcircle of △ACD (blue circle), and then, △ACD is not a Delaunay triangle. Delete △ACD from the undetermined triangles, and △ADE, △ACE, and △CDE are new undetermined triangles.(iii)
PointE is outside and not on the right of the circumcircle of △ABD (red circle), and then, △ABD is still an undetermined triangle.Figure 6
Delaunay triangulation. When pointE is inserted, the undetermined triangles are checked by circumcircle testing. (a) Undetermined triangles. (b) Circumcircle test.
(a)(b)In the DTVC algorithm, instead of storing the Delaunay triangles, we directly calculate the projected volume. The algorithm is shown in Algorithm1. This 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.Algorithm 1: Delaunay triangulation-driven volume calculation.
(i)
Input: the vertex list V(ii)
Output: the volume Vol of V(1)
Vol=0;(2)
Initialize the temporary triangle listT;(3)
Calculate the supertriangle△EFG from V;(4)
V⋅AddE,V⋅AddF,V⋅AddG;(5)
T⋅Add△EFG;(6)
foreachPinVdo(7)
Initialize the edge bufferEdge;(8)
foreach△ABCinTdo(9)
Calculate the triangle circumcircle centerO and radius r;(10)
ifP⋅x>O⋅x+rthen(11)
Vol+= Algorithm 2T,E,F,G;(12)
end(13)
ifOP≤rthen(14)
Edge⋅AddAB,Edge⋅AddBC,Edge⋅AddAC;(15)
T⋅Delete△ABC;(16)
end(17)
end(18)
Delete all doubly specified edges fromEdge;(19)
foreachABinEdgedo(20)
T⋅Add△ABP;(21)
end(22)
end(23)
Vol+= Algorithm 2T,E,F,G;There is a shared algorithm to calculate volume of a set, which is shown in Algorithm1. The calculation of this algorithm is based on equations (12) and (13).Algorithm 2: Triangle projection volume calculation.
(i)
Input: the triangle list T, and the excluded point set P(ii)
Output: the summary volume Vol of T, where the triangle is excluded in which one of vertices is in P(1)
Vol=0;(2)
foreach△ABCinTdo(3)
ifA∉PandB∉PandC∉Pthen(4)
Calculate the projection volumevi of △ABC;(5)
Vol+=vi;(6)
end(7)
end


## 5. Experiments

### 5.1. Datasets

In order to analyze the error of the algorithm, we simulate two point cloud data, as shown in Figure7. Figure 7(a) has a single peak, whereas Figure 7(b) has multiple peaks. The single and multiple peaks have 10.201 and 40.401 points, respectively.Figure 7
Test data for error analysis. The real volumes of (a) and (b) are 2.9284 and 362.6667, respectively. (a) Single peak. (b) Multiple peaks.
(a)(b)To test the volume calculation in real scene, we take 59 photos in a320×110×49 coal shed to generate the point cloud. The resolution of the photos is 3072×2048.


### 5.2. Experiment Setup

In error analysis, we use Matlab R2018b implementation tool. The processor and memory are Intel Core i7-8565U and 16.0 GB, respectively. The operating system is Win10.In the real scene, the 3D reconstruction and volume calculation system run on Ubuntu 18.04 LTS. The CPU, GPU, and memory are Xeon Gold 6226R, Nvidia 16 GB Tesla T4, and 256G DDR4, respectively.The 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:(15)Err=v−vrvr,where v and vr are the measured and the real volume, respectively.


### 5.3. Compare with State-of-the-Art Method

We choose some popular used methods to compare with our proposed method. The notation for these methods is as follows:(i)
Slice-# [29]. Slicing method, and “#” is the slice number.(ii)
P-T3 [17]. The projection method which divides the pentahedron into 3 tetrahedrons.(iii)
DTVC. Our Delaunay triangulation-driven volume calculation.Tables1 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. The results based on the slicing model are difficult to be applied to practical solutions.Table 1
The single peak result.
MethodVolumeError rate (%)Running time (ms)Slice-203.229010.2612.937Slice-503.16818.1913.309Slice-1003.14957.5515.125Slice-10003.07114.8755.491P-T32.92820.0068162.699Our DTVC2.92820.006875.298The real volume is 2.9284.Table 2
The multiple peaks result.
MethodVolumeError rate (%)Running time (ms)Slice-20623.541171.9312.490Slice-50600.361265.5418.017Slice-100596.990664.6124.789Slice-1000543.931349.98172.744P-T3362.66880.00058646.776Our DTVC362.66880.00058286.918The real volume is 362.6667.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.


### 5.4. 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. Figure8 shows a display interface of the system.Figure 8
Volume calculation in a real scene.


## 5.1. Datasets

In order to analyze the error of the algorithm, we simulate two point cloud data, as shown in Figure7. Figure 7(a) has a single peak, whereas Figure 7(b) has multiple peaks. The single and multiple peaks have 10.201 and 40.401 points, respectively.Figure 7
Test data for error analysis. The real volumes of (a) and (b) are 2.9284 and 362.6667, respectively. (a) Single peak. (b) Multiple peaks.
(a)(b)To test the volume calculation in real scene, we take 59 photos in a320×110×49 coal shed to generate the point cloud. The resolution of the photos is 3072×2048.


## 5.2. Experiment Setup

In error analysis, we use Matlab R2018b implementation tool. The processor and memory are Intel Core i7-8565U and 16.0 GB, respectively. The operating system is Win10.In the real scene, the 3D reconstruction and volume calculation system run on Ubuntu 18.04 LTS. The CPU, GPU, and memory are Xeon Gold 6226R, Nvidia 16 GB Tesla T4, and 256G DDR4, respectively.The 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:(15)Err=v−vrvr,where v and vr are the measured and the real volume, respectively.


## 5.3. Compare with State-of-the-Art Method

We choose some popular used methods to compare with our proposed method. The notation for these methods is as follows:(i)
Slice-# [29]. Slicing method, and “#” is the slice number.(ii)
P-T3 [17]. The projection method which divides the pentahedron into 3 tetrahedrons.(iii)
DTVC. Our Delaunay triangulation-driven volume calculation.Tables1 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. The results based on the slicing model are difficult to be applied to practical solutions.Table 1
The single peak result.
MethodVolumeError rate (%)Running time (ms)Slice-203.229010.2612.937Slice-503.16818.1913.309Slice-1003.14957.5515.125Slice-10003.07114.8755.491P-T32.92820.0068162.699Our DTVC2.92820.006875.298The real volume is 2.9284.Table 2
The multiple peaks result.
MethodVolumeError rate (%)Running time (ms)Slice-20623.541171.9312.490Slice-50600.361265.5418.017Slice-100596.990664.6124.789Slice-1000543.931349.98172.744P-T3362.66880.00058646.776Our DTVC362.66880.00058286.918The real volume is 362.6667.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.


## 5.4. 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. Figure8 shows a display interface of the system.Figure 8
Volume calculation in a real scene.


## 6. Conclusion

In the volume measurement of nonconvex point cloud, the triangulation projection method has advantages in accuracy and efficiency. Traditional method of decomposing voxel into tetrahedrons is concise. However, experimental results show that the method of finding analytic solution for voxel volume integral is complex but more efficient. Therefore, the analytic solution method is better than the divide-and-conquer method in efficiency.



---

*Source: 6613264-2021-04-16.xml*