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.
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 Figure 1, 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.
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. 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.
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 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 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 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 Figure 2.
We choose a triangle in the triangular mesh and project it to the base plane to get . The triangles and and the quadrilaterals , , and form a pentahedron, which is called voxel in this paper.
Suppose the voxel set is and the volume of is , and then, the coal volume can be computed as follows:
3.2. Equation of the Top Surface
The equation of the top surface of is required before calculating the volume using integral method. Take the base plane as the plane and z-axis pointing upward to establish the coordinate system. Suppose the top and bottom surfaces of voxel are and , respectively. The coordinates of , , and are , , and , respectively.
We rearrange the order of , , and so that . From the view of the positive direction of z-axis, can be on the right or left side of , 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 , , and .
(a)
(b)
In , , , and are not collinear, so , , and are not in a plane that is perpendicular to the base plane. Thus, .
The equation of the top plane can be expressed aswhich can be transformed as follows:where , , and .
3.3. The Edge Equation of the Projected Triangle
Looking from the positive to negative direction of z-axis, can be on the right or left of line , 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.
(a)
(b)
The coordinates of , , and are , , and , respectively. We first only consider the case of , and then, the equations of , , and can be computed by equations (6), (7), and (8), respectively:where and .where and .where and .
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 can be computed using integral method.
As shown in Figure 4, the projected triangle is divided into two parts by the blue line . It is corresponding that the pentahedron is divided into two parts by a plane that is perpendicular to the plane and passes through the line . The two parts can be calculated by equations (9) and (10), respectively. Then, we can compute the volume of voxel , which is shown in equation (11):
Let , 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 . If or , then there is or , 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 , and is any triangle of . If is a Delaunay triangle of , only the interior of the circumcircle of each triangle in does not contain any point in , 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 in . This section introduces the creation of the supertriangular.
Firstly, find the maximum and minimum coordinates of all points in and directions, which are denoted as , and . These four values construct a rectangle, which contains all the points in . The 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 . According to equation (14), the coordinates of , and are , and , 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.
Thirdly, extend to , which satisfies . Extend to , which satisfies . Then, the coordinates of and are and , respectively.
Fourthly, connect and extend and , which intersect at point . Then, the coordinate of is .
The blue triangle is the supertriangle.
4.3. Triangulation-Driven Algorithm
The key of triangulation is to determine whether a triangle is a Delaunay triangle. In Figure 6(a), , , and are undetermined triangles. When inserting a point , we make circumcircle for each undetermined triangle, as shown in Figure 6(b). According to the relative position of and circumcircle, we can determine the triangle type. The decision rules are as follows:(i)Point is on the right of circumcircle of (green circle), and then, is a Delaunay triangle.(ii)Point is in the circumcircle of (blue circle), and then, is not a Delaunay triangle. Delete from the undetermined triangles, and , , and are new undetermined triangles.(iii)Point is outside and not on the right of the circumcircle of (red circle), and then, is still an undetermined triangle.
(a)
(b)
In the DTVC algorithm, instead of storing the Delaunay triangles, we directly calculate the projected volume. The algorithm is shown in Algorithm 1. 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.
|
There is a shared algorithm to calculate volume of a set, which is shown in Algorithm 1. The calculation of this algorithm is based on equations (12) and (13).
|
5. Experiments
5.1. Datasets
In order to analyze the error of the algorithm, we simulate two point cloud data, as shown in Figure 7. 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.
(a)
(b)
To test the volume calculation in real scene, we take 59 photos in a coal shed to generate the point cloud. The resolution of the photos is .
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:where and 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.
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. The 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.
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. Figure 8 shows a display interface of the system.
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.
Data Availability
The data used to support the findings of this study are available at https://http://www.dropbox.com/sh/mce9q1n9znm8y22/AAB7cCqEPvD4Vl2Z_JxROP_la?dl=0.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Project of Coal Measurement of Huadian Zibo Thermal Power Co., Ltd.