Automatic Stockpile Extraction and Measurement Using 3D Point Cloud and Multi-Scale Directional Curvature

: With the rapid increase of power supply demand, a large amount of stockpiles of coal have been formed during the process of coal excavation and transportation between the mines, ports, power plants and etc. Quantitative parameters, especially the volume of stockpile are important for the planning of coal production and consumption. Although laser scanning can collect the dense 3D coordinates of the stockpile surface for its quantiﬁcation, stockpiles of coal have irregular shape, size, height, and base conditions, and may be overlapped with each other, which makes it hard to measure the different stockpiles automatically and accurately. This paper proposes an algorithm to extract and measure the stockpiles from the 3D point cloud data using the multi-scale directional curvature. Firstly, the second-order directional curvature analysis along multiple directions and at multiple scales is proposed to extract the distinctive ridge of crest in the point cloud of stockpiles. Then, starting with the crest points, a competitive growing strategy is proposed to accurately locate the points of slope in the stockpile. Finally, the stockpile’s volume is calculated by reconstructing the complete points of crest and slope with many meshes and triangular prisms through the subsurface ﬁtting and surface reconstruction. Experiments on both the synthetic and real point cloud of stockpiles demonstrate that the proposed algorithm can extract the stockpiles with the average accuracy over 93.5% and measure the volume of stockpiles with the average precision over 93.7%. It is promising for automatically measuring the stockpiles like sand, soybean, etc., and facilitating the scientiﬁc storage and transportation management, as well as the maintenance of safety inventory during operation.


Introduction
As power supply demand grows rapidly, the demand for coal is also increasing. In the production and transshipment process of coal, stockpile is the basic unit for management. A large amount of stockpiles of coal have been formed while they are excavated and shipped at the mines, ports, power plants, etc. Since stockpiles take huge space, it becomes hard and time-consuming to monitor/measure them accordingly, which brings difficulty to the planning of coal production and consumption. It is necessary to quantitatively measure the coal stockpiles to facilitate the scientific storage and transportation management, as well as the maintenance of safety inventory during operation. This paper is then motivated to automatically reconstruct the stockpiles and measure their volume.
Commonly, the volume of the stockpile can be calculated once loaded on trucks or trains or from their weight. Besides, terrestrial geodetic methods using total station are often used to obtain information about the stockpile, which need to walk over the stockpile to manually sample many control points on the stockpile to estimate the volume. However, the above methods are time-consuming and laborious, and the precision is low. More effective and reliable methods for stockpile measurement should be utilized.
With the development of sensors and computer technology, three-dimensional (3D) point cloud data was involved to achieve better measuring accuracy and more detailed spatial morphology of the stockpiles. 3D point cloud collect the dense 3D coordinates of the stockpile surface, quantitative parameters of stockpile such as height, floor area and volume can then be computed through the processing of 3D point cloud data.
There are many ways to obtain 3D point cloud data, such as depth camera, stereo camera, laser scanner, etc. In some studies [1,2], depth cameras based on structure light or Time-of-Flight were used for 3D reconstruction. Combined with many denoising approach, the accuracy of 3D reconstruction and 3D point cloud data obtained by depth camera in real-time is acceptable. it is an easy-to-use and low-cost way to use depth camera, however, the effective distance of depth camera tend to be so short that depth camera is unsuitable to be performed in large field [3] in the coal yard, and the sunlight can also impact the precision of depth camera under outdoor conditions [4]. A method for 3D reconstruction using stereo vision camera pair has been described in [5], by which objects can be reconstructed from the dense matching. The work in [6] combined the efficient stereo matching with a multi-view linking scheme for generating consistent 3d point clouds, allowing real-time 3d reconstructions from large-scale imagery. The work in [7] used UAV with stereoscopic aerial photographic pairs for capturing point cloud data and estimating stockpile volumes. 3D point cloud obtained by stereo vision technique have good quality and high resolution. But when it comes to some scene where the image features are not obvious, such as ground, coal, desert, etc., stereo vision technique is hard to work, as image matching, which is the most critical step for stereo vision, can't be achieved effectively if there exist little distinctive features.
Compared with image-based photogrammetry methods, or the manual survey using GPS and total station, laser scanning technique has better accuracy, robustness and efficiency. It is extensively used for topographic measurements in recent years. The work in [8] estimated the errors and accuracy of laser data acquired by various laser scanning systems for topographic applications and the results showed that laser scanning can be used for the extraction of topographic information for various purposes. In [9,10], 3d Laser Scanning and GPS technology were used to acquire landslide point cloud data and to compute earthwork volume. In this paper, we take advantage of 3D point cloud data from the laser scanning to compute key parameters of stockpiles.
To measure such parameters as volume, height of stockpiles from the 3D point cloud data, work [7,11] studied cases in an open pot quarry and stockpile using Pix4D software. It achieved high accuracy of measurement and cost much less time than manual surveys. The work in [12] also calculated the volume of stockpile using 123d software on point cloud data, and resulted in a higher speed and accuracy of volume measurement than conventional methods.
Although volume calculation methods using point cloud have better speed and accuracy than conventional methods, there are still challenges to obtain accurate quantitative parameters using 3D point cloud data. As the basic unit for the management of the actual production and transshipment process, different stockpiles may represent different production date, method, type, or transportation mode, and need to be distinguished accurately on the boundary even if they are overlapped between each other. Most of the existing researches need to delineate the boundary manually on the point cloud. Stockpiles of coal have irregular shape, size, and base conditions, and may be overlapped with each other, which make it hard to identify the boundaries clearly and measure stockpiles separately. Therefore, there is a need to automatically extract and measure the stockpiles from the 3D point cloud data.
Extracting the stockpiles from the whole scene can be regarded as a point cloud segmentation problem. Currently, the existing point cloud segmentation algorithms include region growing algorithms, model fitting algorithms, edge-based algorithms, etc.
(1) Region growing algorithms starts from seed points and then classifies the points with similar characteristics into one class, such as curvature [13,14], point normal direction [15], local convexity and dimension features [16], etc. And the work in [17,18] used an octree-based region growing method according to a coarse-to-fine concept, which improved the segmentation precision. However, these methods are sensitive to the location of initial seed points, and the inaccurate selection of seed points will affect the segmentation results. Moreover, these algorithms rarely consider the overgrowing between overlapped stockpiles.
(2) Model fitting algorithms can also segment point cloud accurately. In work [19,20], an automated detection of 3D geometric shapes from point cloud data was realized based on Hough Transform. The work in [21] presented an iterative sphere fitting algorithm, which can extract both single spheres and multiple spheres. The work in [22] completed cylindrical objects estimation in point cloud. It used RANSAC to estimate the center axis and radius using the trace of a sphere model, and confirmed of the constructed center axis data based on PCA. However, model fitting is mainly efficient for the detection of geometrically simple and regular shape such as cylinders, spheres, cones, cubes, etc. For stockpiles of coal, their shapes may be irregular, or even two of them may be overlapped, which is difficult to describe with any specific geometric models, so that the model fitting of the entire stockpile is not suitable for extracting it.
(3) Edge-based methods firstly detect edges to outline the borders of different regions according to point cloud features, such as normal vector orientation, elevation difference, etc., and then group the points based on the boundaries. The work in [23] realized road boundaries extraction through detecting the curbs in real-world urban scenes, where there is either an elevation jump or a gradual elevation change. The work in [24] put forward a hybrid method based on region growing and model fitting for edge detection by analyzing angular gap of neighborhoods. However, those method can only preform when the edges demonstrate the abrupt changes between the normal directions or elevations of adjacent surfaces, which is not suitable for detecting the continuous stockpiles of coal.
As a matter of fact, although a stockpile has three parts of crest, slope, and base, it is distinctive and perceptible for the crest, which can then act as the starting point to extract and measure the stockpile. However, due to the various height, shape, size, and base conditions, it is hard to locate the crest of different stockpiles. This paper first proposes an innovative perspective of directional curvature to extract the distinctive ridge of crest in the point cloud data. It characterizes the ridge using the curvature analysis along multiple directions and at multiple scales, resulting in invariant to the height, shape, size, and base conditions of different stockpiles. Then, starting with the crest points, a competitive growing strategy is proposed to accurately locate the points of slope in the stockpile, no matter whether there exist flatland, overlapping on the slope or not. Finally, the stockpile's volume is calculated by reconstructing the complete points of crest and slope with many meshes and triangular prisms, which is simple yet general enough for different stockpiles.
The main advantages of this paper are as follow: 1. The proposed algorithm has great robustness not only to different stockpiles of various height, shape, and size, but also to different base conditions of various terrains where the stockpiles are located, as it can adaptively extract the crest and boundary of the stockpiles with the proposed directional curvature. 2. The proposed algorithm can achieve precise result of volume measurement for different stockpiles, by representing the structural features of stockpile into the volume element of triangular prisms. 3. The proposed algorithm can identify the overlapping between stockpiles with the competitive growing strategy, resulting in effective quantification of overlapped stockpiles.

Methods
As shown in Figure 1, the proposed method mainly includes the following steps: 1. Crest Extraction: Multi-scale 3D directional curvature is proposed to locate the crest points among the point cloud of stockpiles, with the help of the preprocessing and voxelization of the point cloud. 2. Stockpile Extraction: Based on the gradient calculation, the competitive growing strategy is developed to extract the slope points and delineate the boundary of stockpiles, resulting in the accurate and complete points of the stockpiles. 3. Volume Calculation: The stockpile's volume is precisely calculated by reconstructing its point cloud with many meshes and triangular prisms that are obtained through the subsurface fitting and surface reconstruction.

Registration
For a large yard of stockpiles, it is hard to collect its complete point cloud by scanning only once. A stockpile could be occluded by another one or other facilities in the yard if the scanner is not installed highly enough. The back-side slope of a stockpile is usually invisible to the scanner due to the occlusion of the front-side slope. Therefore, multiple scans are employed in practice to partially collect the point cloud of a stockpile at different positions and aspects. As shown in Figure 2, there exists overlapping between the two neighboring scans. And some distinctive targets like the crane in Figure 2 can be used to calibrate the transformation between neighboring scans, so that the partial point clouds can be registered together to obtain the overall surface of a stockpile.  For the two neighboring scans in Figure 2, where P and Q are the coordinates of all the points in two scans, the correspondent target points q i ∈ Q in Figure 2b q i are the center of points p i and q i respectively. To stitch the two scans together by registration, the relative coordinates of the target points with respect to the center are first obtained as follows: where N is the total number of corresponding target points in the overlapped area. p i records the relative translation between p i and the center of target points in a scan, and ensures the diverse distribution of them. The correct rotation parameters between two scans should achieve the minimal error of reprojection, if they are expressed together within one common scanning system, i.e., where P = [p 1 , p 2 , . . . , p N ] and Q = [q 1 , q 2 , . . . , q N ], and P Q T can be decomposed with the single value decomposition (SVD), i.e., where V and W are orthogonal matrices, and S is a diagonal 3 × 3 matrix. From the property of Tr(AB) = Tr(BA), we obtain that Tr(R r P Q T ) = Tr(SW T R r V). Because the absolute value of each element of the orthogonal matrix W T R r V is no more than one, Tr(R r P Q T ) achieves the maxima if W T R r V = I holds. Since V and W are orthonormal, the rotation matrix R r between two scans can then be calculated as And the translation between two scans equals to the center translation of target points, i.e., Based on the rotation and translation from Equations (4) and (5), the point cloud P of the first scan can be stitched together with Q, and updates the overall point cloud P as follows, As shown in Figure 3, where the blue-and red-colored point cloud denote the two scans of P and Q respectively, both the crane target and stockpile become well registered in the stitched update of P, and the whole yard/stockpile appears incrementally by stitching the scans one by one.

Normalization
Due to the variety of incident angle while scanning, the stockpiles are possibly not erected uprightly with respect to the ground in the registered point cloud. As a result, the Z-axis in Figure 3 does not align with the elevation increase of the stockpile, which makes the volume calculation complex. Therefore, the registered point cloud of the yard is normalized to align the ground with the XOY horizontal plane.
(a) before normalization (b) after normalization Firstly, with the ground points in the yard, we can get the normal n z of the ground plane by 3D plane fitting. The two orthogonal axis n x and n y can then be chosen on the ground plane as Figure 4a, so that n x and n y are parallel to the boundary of the yard. The rotation matrix R n for the normalization can be obtained as follows: which describes the rotation between the two coordinate systems o − xyz and o − x y z in Figure 4a,b respectively. After normalizing with R n , the point cloud P of stockpiles is standardized as follows: The Figure 5 shows the same stockpiles as Figure 3, where the point cloud is colored according to the magnitude of z-value. As shown in Figure 5, stockpiles in the registered point cloud of Figure 3 becomes much perceptible above the ground after the normalization.

Crest Extraction
After the preprocessing, a stockpile looks like Figure 6a, which can be divided into three parts: the crest at the top, the bottom base,and the transitional slope. They are colored in red, green, and blue respectively. Some flatland may appear in the slope if the stockpile becomes large enough upon the base. A stockpile is distinctive and perceptible for the crest, which can then act as the starting point to extract and measure the stockpile. However, due to the various height and shape, it is hard to set a simple elevation threshold to locate the crest of different stockpiles. Although the gradient at the crest ridge approximates zero, the gradient at the base is also close to zero and becomes indistinguishable with the crest. Indeed, the crest point A is distinguished from the slope point B and base point C, in that the curvature is large and negative at the crest ridge of the stockpile, as shown in Figure 6b,c.
In contrast, although the slope point B has a large gradient, its curvature is close to zero. There exists sudden change of shape at the base point C, but it forms a valley and has a positive curvature. This paper is motivated to take advantage of the curvature to extract the stockpile crest first in this section, followed by the growing of its slope in the next section.

Voxelization
To calculate the curvature, the dense point cloud of the yard is divided into regular grids on the normalized XOY plane. Figure 7a shows the original point cloud. The Figure 7b is the result of voxelization of Figure 7a by dividing the yard point cloud into grids on the XOY plane, where grid size S is fixed and identical for the X and Y axis. All points located within a grid are represented with their central point in the 3D space, so that the 3D point cloud are transformed into 2D elevation image. And for the grid with no laser point, the elevation is inferred with the bilinear interpolation of the neighboring pixels on the XOY plane of the elevation image. As shown in Figure 7b, each grid is colored according to the elevation. The redder a grid is, the higher elevation it is located at. With the elevation image, both the gradient and curvature at a space point can be readily calculated due to the regular organization of 3D points.

Multi-Scale 3D Directional Curvature
As shown in Figure 8, the curvature of a space point varies with the direction. The black line in the first column of Figure 8 represents the direction along which the stockpile is sliced and the curvature is calculated. It can be seen from the second column that the directional profile of the stockpile is quite different with the varying slope along different directions of the stockpile. Near the crest point, since the vertical slope is much flatter than the horizontal one, the stockpile's profile changes more slowly along the vertical direction than the horizontal direction, resulting in a larger curvature at point A along the vertical direction. The images of derived multi-scale curvature along two directions are shown in the last column of Figure 8. It can be seen that the complete crest can be extracted by integrating the results along different directions, although the horizontal curvature of point A in Figure 8d is not as strong as its vertical curvature in Figure 8h. Therefore, for a pixel in the elevation image, its curvature is calculated along multiple directions for the crest extraction. Assume that the elevation image is represented as z = f (P(x, y)) with z representing the elevation of the grid P(x, y), the directional curvature K at P(x, y) along direction θ can be obtained as follows: where ∂z ∂θ | P(x,y) and ∂ 2 z ∂θ 2 | P(x,y) represent the first and second order directional derivative of elevation z at P along direction θ.
To calculate the curvature more accurately for various changes of the profiles in the second column of Figure 8, the multiple Gaussian scale spaces are employed to adapt to the local curvature of the profiles. As shown in the third column of Figure 8, compared with the rapid change of crest profile in the first row, a larger scale can be chosen for the second row while calculating its crest curvature due to the slow change of profile. To adaptively calculate the crest curvature of stockpiles of various slope size and shape, the directional curvature is then formulated within the multiple Gaussian scale space as follows: where G (σ,θ) is Gaussian smoothing filter. The σ is used to adapt to the scale, , n = 1, 2, · · · , N, a total of N scales are analyzed for the curvature of a point. The directional derivative G filter at the direction θ is synthesized by taking a linear combination of G σ,0 • and G σ,90 • : Where G 0 • = ∂G ∂x , G 90 • = ∂G ∂y . And a G filter at an arbitrary direction θ is synthesized by taking a linear combination of G 0 • , G 60 • and G 120 • using an interpolation function k(θ, θ j ) [25]: where . With the multi-scale 3D directional curvature analysis of Equation (10), the curvature of a space point is analyzed along multiple directions and with multiple scales. The final curvature K P of the point P is decided by the direction and scale with which the maximal absolute curvature is obtained among all the directions and scales, i.e., where σ * , θ * = argmax Therefore, the point becomes the candidate of crest point provided that its final curvature is negative and its absolute curvature is big enough, i.e., where C is negative and denotes the curvature of the minimal perceptible stockpile.
Since the crest points are located at a continuous area in the yard, the morphological closing operation could be added to fill the discrete hole of noise and refine the crest area after the multi-scale 3D directional curvature analysis.

Slope Extraction
Besides the crest, the slope also accounts for a large portion of a complete stockpile. It is spatially connected with the crest. With the result of crest extraction, the slope of a stockpile can then be located by growing the crest over the 3D stockpile surface. The slope is differentiated from the base because it is steeper and its gradient is much larger. Therefore, the growing process is based on the gradient and stops at the base or the overlapped areas of multiple stockpiles.

Calculation of Spatial Gradient
For 3D points located in the grid P i of the elevation image, their gradients can be assumed identical and equal to the degree of elevation change. As shown in Figure 9, there are 8 different gradients along a grid's neighboring directions, such as J, J ∈ 1..8, which are defined as follows: where (x i , y i , z i ) and (x j , y j , z j ) are the 3D spatial coordinates of P i and its neighbor along the direction J, and G (P i ,J) reflects the elevation change of P i along the direction J.
For stockpiles of different shapes, their gradients vary along different directions. The cone-shape stockpile B cannot be perceived in Figure 10a by the gradient direction of J = 4, but becomes obvious with the gradient direction of J = 5 in Figure 10b. It is because the stockpile B is not strictly symmetric and drops faster along direction J = 5. For the ridge-shaped green stockpile C, it spreads anisotropically, thus its gradients along direction J = 4 and J = 5 are both obvious in Figure 10. As a result, all the 8 gradients of a point are considered for judging whether it is located on the slope.

Slope Growing
Starting from the crest, the growing process of slope could be like either Figure 11a or b, which respectively represents two kinds of slope: without or with flatland inside the slope area. The growing continues as long as the gradient keeps large enough, and stops at the base. The slope grows along all the spatially neighboring directions so that it can adapt to the different shape of stockpiles. Even for the flatland inside the slope, the growing can still move around and avoid the early stopping.
Initially a queue of crest points is prepared as the seeds of slope growing. The immediate neighbors of the seed points are visited clockwise and checked if they are steep away from the seed. If the directional gradient of a neighbor is larger than a threshold, it is marked as slope and will serve as a new seed by pushing into the queue of seeds for the next growing; otherwise, the directional gradient of the next neighbor will be checked until all the neighbors of the seed have been checked. The seed is then removed from the queue once all its neighbors have been checked. The gradient check process continues so long as the queue is not empty. Finally the growing stops at the absolute flat base area in Figure 11a  For the flatland inside the slope in the left of Figure 12, although the growing stops along the red arrow due to the flatland, it moves around the flatland along the green arrow and finally reaches the blue base. It is because the gradient along the green arrow is still large enough to keep the growing around the flatland. As a matter of fact, though there is flatland inside the slope, the stockpile is continuous in the 3D space and will reach the ground of the yard. The flatland is always surrounded by areas of high gradient of slope that will grow towards the ground. Therefore, the flatland creates a closed not open hole for a stockpile, and can be filled with the morphological closing of holes on the XOY plane of the elevation image, as shown in the right of Figure 12, where the gray flat area surrounded by the green grown area is filled, and the complete point cloud of the stockpile is obtained.

Boundary Delineation
During the actual production and transshipment process, stockpile is the basic unit for management. Different stockpiles may represent different production date, method, type, or transportation mode. Therefore, for scientific storage and transportation, we distinguish clearly between different stockpiles, although the volume of stockpiles won't change if they are not separated from each other.
As shown by Figure 13a, there appears over-growing of the red stockpile into the blue stockpile, which is caused by the overlapping of the two stockpiles and the unavailability of base in this case. Basically, the growing process proceeds faster along the direction of larger gradient, thus the overlapped green stripes area are firstly visited by the red stockpile because it is much steeper than the blue one. Actually some points located in the green diagonal striped area are much closer to the blue stockpile and should be attributed to the blue one. Stockpile often accumulates from bottom to top. In practice, separate stockpiles are close to symmetrical cones, and the centers of stockpiles are unbiased. A competitive strategy of growing is adopted for the overlapped points, which follows the principle of "the closest the first", i.e., the overlapped point is assigned to stockpile whose center is much closer to the point. As shown in Figure 14, the center of the stockpile is updated from the original K 1 to the new K 2 after one iteration of slope growing. It is dynamically calculated by the average coordinates of all points belonging to the stockpile. During the process of competitive growing, the minimal distance to all the stockpile's center is therefore recorded for each point besides its membership to the closest stockpile.  As a result, some overlapped points will be re-assigned to blue stockpile if they are much closer to the crest of blue stockpile, although they are previously assigned to the red one, as shown in Figure 13b.

Volume Calculation
After extracting the crest and slope, the complete point cloud of a stockpile is ready for the volume calculation. But the stockpile has arbitrary shape and follows no regular model of calculation. Although the previous voxels are ready for the volume calculation, they do not match the smooth and continuous surface of stockpiles well. Therefore, the triangular meshes are constructed over the voxels to approximate the surface as much as possible. After reconstructing the stockpile with many meshes, the stockpile's volume is calculated by summing up the volume of triangular prisms that are formed between the meshes and the base. First, the base plane around the stockpile is obtained by fitting the boundary of stockpile into a subsurface. Then, the manifold surface of the stockpile is reconstructed with many triangular meshes. Finally, the entire entity of stockpile is decomposed into triangular prisms by projecting the meshes onto the base plane. Adding the volume of triangular prisms together results in the volume of stockpile.

Subsurface Fitting
In this paper, the subsurface is fitted based on the eigenvalue method using the boundary points of the stockpile. The overlapped boundary points of two stockpiles will be removed in the fitting, because overlapped points is generally not at the intersection of stockpiles and subsurface, and doesn't represent the base of the stockpile. Let a spatial plane be expressed as follows: where a, b, c is the unit normal vector of plane, i.e., a 2 + b 2 + c 2 = 1, and d is the distance from the origin of coordinate system to the plane. Assume there are n boundary points (x i , y i , z i ), i = 1, 2, . . . , n for the plane fitting, based on Equation (16), the distance from a boundary point (x i , y i , z i ) to the base plane can be calculated as follows: The best fitted base plane should minimize the total distance ∑ i d 2 i between all boundary points and the base, subject to the condition of a 2 + b 2 + c 2 = 1. Solving the conditional extremum problem with Lagrange's multipliers, we have  (18), the problem of solving a, b, c can be transformed into calculating the eigenvalue and eigenvector of the equation Ax = λx. Therefore, the eigenvector corresponding to the minimum eigenvalue of A is the coefficients a, b, c of the base plane. And the distance d of the base plane is obtained as follows:

Surface Reconstruction
To approximate the smooth surface of stockpile as much as possible, we construct a triangular mesh of the stockpile surface with the voxels, as shown in Figure 15, where the underlying grids are obtained from the voxelization. But the stockpile has arbitrary shape and follows no regular calculation. For the irregular shape of stockpiles, their point cloud are decomposed into many regular triangular prisms for the volume calculation. Based on the previous voxelization of point cloud, we further divide the four adjacent points of a quadrangle on the mesh into two triangles. In order to make a more accurate measurement of the volume, the quadrangle is always divided along the diagonal points with small difference of elevation. As shown in the right of Figure 15, the quadrangle mesh of ABCD is divided into two triangular prisms along the ACC A , because the elevation difference between A and C is smaller than that between B and D, i.e., |Z a − Z c | < |Z b − Z d |. In this division, the resultant prisms become much approximate to the stockpile's surface spanning by the point cloud.

Volume Calculation
After the surface reconstruction, stockpile is divided into many regular triangular prisms. The stockpile's volume can then be calculated by summing up the volume of all triangular prisms. Each triangular prism, as shown in Figure 16, is composed of the surface points A, B, C and their projection on the fitted subsurface base A , B , C . The volume of the triangular prism can be calculated as follows: where S A B C is the area of triangle A B C on the base plane. The total volume of the stockpile can be obtained by summing up the volume of all the triangular prisms.

Experimental Results
In order to evaluate the performance of the proposed algorithm for extracting and measuring stockpile, a diverse set of experiments are given with point cloud from both the synthetic and real yards of stockpiles. Stockpiles in the real yard could be very large and have arbitrary shape. It is hard and time-consuming to validate the algorithm accuracy for various stockpiles in the real yard. Therefore, a simulation device was designed which is able to collect point clouds of stockpiles of different shapes and sizes, and measure the ground truth of volume for the stockpiles. Experimental results with the synthetic and real data are evaluated both qualitatively and quantitatively to demonstrate the performance of the proposed algorithm.
The proposed algorithm is implemented in C++ language, and results of each algorithm step are also rendered carefully for better qualitative evaluation. For the voxelization of the point cloud from real yards, the grid size S is set according to the point density so that we can adaptively achieve a high precision of volume measurement without compromising the computation time cost. A large voxel size will speedup the computation but reduce the precision of volume, because the surface of stockpile becomes roughly sampled in this case. With a small voxel size, a high Precision of volume measurement will be achieved but it is time-consuming. Therefore, the choice of S is a balance between the precision and time, and keeps fixed to include around 20 points within a grid for real experiments. The curvature of a point is analyzed with 7 scales of 7,9,13,19,25,35, 49 pixels and 4 directions of 0, 45, 90, 135 degree. The curvature threshold C is set according to the minimal perceptible stockpile (e.g., the stockpile of a 30-degree slope at the minimal scale) in the yard, and the gradient threshold of slope growing equals to 0.1 so that the growing will stop at the minimal slope of around 5 degree.

Performance Measure
To measure stockpiles that are of arbitrary shape and size, and located at varying subsurface base,the proposed algorithm includes three steps: crest extraction, slope extraction, and volume calculation, whose results are evaluated separately. For the extracted crest, whose number is small in most yards, it is easy to check whether there exist false positives and negatives. Thus, only qualitative rendering of crest point clouds is given for growing the slope of stockpiles, since the number of them has nothing with the growing. But the number of the complete point cloud of a stockpile is decided by the performance of extraction, and affects the precision of volume measurement. Therefore, the proposed algorithm is evaluated quantitatively in terms of accuracy of stockpile extraction and precision of volume calculation.

Accuracy of Stockpile Extraction
The accuracy of stockpile extraction from the point cloud is important for the stockpile measurement. To evaluate the extraction accuracy, the intersection-over-union (IoU) is employed as follows: where A and A are the set of points that are extracted by professional operators and the algorithm respectively. IoU measures the degree to which the algorithm mimics the operator in terms of extracting point cloud of stockpile. The higher IoU is, the more accurate the extraction result coincides with the ground truth. The ground truth for the evaluation is carefully set up manually with the point cloud selection tool of ArcGIS 10.2.

Precision of Volume Calculation
The precision of volume calculation is measured by the relative error of the algorithm compared to ground-truth volume. For the synthetic data, the ground truth is obtained with the volume measuring cup in ahead of creating a stockpile. The volume of stockpile in real yard is also carefully obtained by transportation vehicles with volume metering capability.
where V and V denote the ground-truth and calculated volume of the proposed algorithm respectively for a stockpile. The higher Precision-of-Volume (PoV) is, the more precisely the algorithm measures the volume of a stockpile.

Experiments with Synthetic Data
To comprehensively evaluate the proposed algorithm, a diverse sets of synthetic point cloud with known ground truth are created with a simulation device. They include one or more stockpiles, and/or are overlapped with each other, and/or located at an uneven subsurface base in the yard. Certainly, the shape and size of the stockpile vary a lot with every creation.
As shown in Figure 17, the simulation device for collecting point clouds of arbitrary stockpile is inspired by [26], which is composed of three parts: the Kinect V2 depth sensor, a projector and sandbox of stockpile(s). The whole device is placed on a horizontal level. The Kinect is installed downward and 2 m above the sandbox. It captures the distance of sandbox away from its imaging center, resulting in point clouds of the stockpiles in the sandbox, with the average depth resolution around 2 mm. The Kinect V2 has a depth image resolution of 512 × 424 pixels with a field of view of 70.6 × 60 degrees. Because of this installation, the complete stockpiles' point clouds can be collected with the near and downward imaging. There is no need of stitching point cloud of multiple scans for a stockpile and then voxelizing point cloud. The projector is used to vividly render the stockpile by its volume, elevation, etc. once it is calibrated well with the Kinect depth sensor. Since the sand can be shaped arbitrarily, it is helpful to create stockpiles of various shape, size, overlapping, and subsurface base within the sandbox. For each stockpile, its volume can be readily metered in ahead of the creation, through which we have the ground truth regarding the stockpile's volume. Once the point cloud of the stockpile is collected by the Kinect, the proposed algorithm is called to extract and measure it automatically. The performance of the algorithm is quantitatively evaluated with the ground truth.

Experimental Results with Typical Yard
As shown in Figure 17, the sand that fills up the 5 litres measuring cup is put into the sandbox and under the Kinect depth sensor. The point cloud of the sandbox looks like Figure 18a, where the ridge of stockpile is obvious due to its protruding elevation. It is the typical yard with one significant stockpile in the middle.
After we get the stockpile's point cloud with the Kinect, the crest and slope extraction results are shown in Figure 18b, where the reddest points in the middle are the crest points and the red points located between the gray base and the crest denote the slope points. It can be seen that both the crest and the complete stockpile are extracted well. Figure 18c is the top view of Figure 18b to better demonstrate the accuracy of boundary delineation of the proposed algorithm, where the black edge between the slope and the base is delineated by a professional operator. It can be seen that the boundary of the slope growing matches well with the professional delineation. The little difference in the rectangle A and B of Figure 18c could slightly reduce the IoU measure of slope extraction. But since they are located on the boundary, the elevation of A and B is very close to the base and cannot form perceptible prism at all, contributing little to the stockpile's volume measurement.  The quantitative result of stockpile extraction and measurement of typical yards is given in Table 1 for the proposed algorithm, where 10 experiments of stockpiles of different sizes are tested, and the IoU and precision are calculated for each experiment in the second and last column . From Table 1, we can see that the average IoU and precision are 93.9% and 96.2%. The minimal extraction accuracy of stockpile point cloud is over 90%, and shows good matching with the human delineation. Due to the limit of Kinect's resolution of depth measurement, the precision of volume measurement fluctuates a little, which relatively brings more impact on the extraction and measurement of small stockpile than on large stockpile. The minimal measurement precision is over 85%. It demonstrates that the proposed algorithm can measure stockpiles of different sizes precisely. Table 1. Performance of stockpile extraction and measurement in typical yard.

Accuracy of Extraction Precision of Measurement
IoU/% Ground-Truth Volume/L Calculated Volume/L PoV/%

Experimental Results with Diverse Yard
Besides the typical stockpiles, there is possibly flatland in the slope, or uneven subsurface surrounding the stockpile, which brings difficulty for the crest or slope extraction by elevation. In addition, the occurrence of multiple stockpiles in a yard, possibly overlapped with each other, requires accurate location of the stockpile boundary at the end of slope growing. Therefore, a diverse set of yards are investigated to show the performance of the proposed algorithm in this section (1) stockpile with flatland slope As shown in Figure 19a, there exists flatland in the rectangular area of the slope, where the gradient becomes close to zero although it is located on the slope. The gradient-based slope growing could stop early at the flatland instead of the subsurface base of the stockpile.
The final result of slope growing is shown in Figure 19b, where the flatland doesn't stop the location of the stockpile boundary and is successfully filled up. Since the slope points can be grown along all directions of the flatland boundary as long as the gradient is above the threshold, the flatland only stops partial directions of growing through it and enables the growing around it. Therefore, with the proposed omnidirectional spatial growing of slope, the slope points are completely extracted even there are flatland inside, as shown in the top view of flatland in Figure 19c. (2) stockpile on a slant When the stockpile is located on a slant, it is difficult to extract it by elevation, as shown in Figure 20a, where the points are rendered by the elevation. The result of crest and slope extraction is shown in Figure 20b. Both the crest and slope are accurately extracted with the IoU of 95.9% and the precision of 97.8% for the volume measurement. Actually, after the normalization of point cloud, the effect of uneven base on the stockpile extraction is alleviated a lot. Additionally, the proposed multi-scale directional curvature analysis works at the change of slope gradient, thus it is robust to slant base during the crest extraction.   Table 2, which achieves a high IoU of 94.4% and precision of 96.5%, demonstrating that the proposed algorithm performs well with various shape, size, and location of the stockpile.  (4) multiple stockpiles with overlapping between each other For the first column of Figure 22, there are overlapping between different stockpiles. In this case, the gradient-based slope growing could reach the overlapped boundary at different speeds, causing over-growing of the boundary of one stockpile into the other one's. With the proposed strategy of competitive growing, points located at the overlapped area are assigned to the crest that is much closer to them. The extraction results of two, three, and four overlapped stockpiles are shown in the second column of Figure 22, where no over-growing appears between stockpiles. It can also be seen from the third-column top view of the boundary in Figure 22 that the intersecting boundary of differently-colored stockpiles matches well with the professional delineation.   Table 3 gives the quantitative results of Figure 22, from which we can see that the proposed algorithm performs well for stockpiles of overlapping. The average IoU of extracting complete points of stockpile achieves 93.5%, and the volume is well measured with a precision of 93.7%. Table 3. The quantitative results in overlapping stockpiles.

Terrain Condition
No.

Experiments with Real Data
Besides the synthetic experiments with the sandbox, real experimental results at a coal yard are given in this section. The size of the yard is 14904 m 2 . There are four stockpiles of different shape and size inside the yard. The point cloud of the stockpiles and yard was collected by the SICK LD-LRS3611 which can work 250 m away the outdoor yard with the angle resolution of 0.125 degree and range resolution of 3.8 cm. Tested in the paper [27], it was proved that this laser can obtain a very accurate range data of coal stockpiles even in the environment with dust, and is suitable for practical use in our experiments. To get the complete yard, it was collected using 15 scans at different positions/orientations and stitched together as a whole.
The ground truth boundary of the stockpile is then elaborately delineated by professional operators with the stitched point cloud of the yard, and the ground truth volume of the stockpile is obtained from the weighting outcome of professional measuring vehicle before the coal transportation.

Results of Crest Extraction
After the preprocessing, the multi-scale directional curvature is calculated in Figure 24a. It can be seen that there are four blue areas with small curvatures after the the multi-scale directional curvature analysis. They exactly correspond to the crest of four stockpiles. To investigate the effect of the proposed curvature analysis, the results of the red rectangular area in Figure 24a are detailed as an example in the Figure 24b-e, where the first and second column are the result of vertical and horizontal curvature analysis; the first and second row denote the curvature at a scale of 49 pixels and 7 pixels, i.e., 34.3 m and 4.9 m respectively. The red rectangular area in Figure 24a is located at a large stockpile spreading horizontally in the local neighborhood. It can be seen that the crest of the rectangle are not well located with the small scale of 7 pixels in Figure 24d,e, but become perceptible with the large scale of 49 pixels in Figure 24b,c. Meanwhile, the crest of the rectangle in Figure 24b becomes more obvious with the vertical curvature analysis than the horizontal one in Figure 24c. Therefore, the curvature of the rectangular area is decided by the vertical curvature at a large scale as Figure 24b. As long as the curvature is greater than a threshold, the point is regarded as a crest point.
In the experiments, the curvature threshold C is set to −0.035 according to the minimal perceptible stockpile of 7 pixels, and the result of crest extraction is shown in Figure 25, where the curvature of a point is analyzed with 7 scales and 4 directions. It can be seen from the colored points that the crests of four stockpiles are well located after the multi-scale directional curvature analysis.

Results of Slope Extraction
Based on the result of crest extraction, the slope of the stockpiles are grown in Figure 26. The gradient threshold of the slope growing is set as 0.1, meaning that the growing will stop at the boundary of a slope of less than 5 degree. From the differently-colored points among the gray base of the coal yard, the complete points of the stockpiles are extracted well. The boundary of stockpile a and b is divided without over-growing due to the employment of competitive strategy during slope growing, as shown in Figure 26b. The quantitative results of the stockpile extraction and measurement are listed in Table 4. For the stockpiles in the real yard, the average IoU of stockpile extraction is 93.9%, and precision of volume measurement achieves 94.6% on average. It demonstrates that the proposed algorithm performs well for stockpiles of different size and shape in the real yard.

Comparison with Other Methods
Compared with the large yard, the stockpiles occupy a little. In order to segment the stockpiles, the intuitive method is based on the underlying ground of the stockpiles, which can be obtained with the median filtering. The window size of median filtering is set to 1.5 times the maximum stockpile. The clusters of points whose height is larger than the ground are regarded as stockpiles. As shown in Figure 27, the stockpile a, b, and d in the real experiments are partially extracted while the stockpile c is missed. The average IoU is 55.9% for the stockpile a, b, and d, which is smaller than the 94.2% by the proposed algorithm. Compared with Figure 26, it is hard to choose a right window size that is suitable for all the stockpiles of different sizes in Figure 27. As a result, the underlying ground of the stockpiles could be over/under-estimated. Also, it is hard to obtain the underlying ground if there exists overlapping between stockpiles. Instead of the bottom-up location of stockpiles, the proposed curvature analysis locates the points of stockpiles from up to down, which is more robust to the different shape, size, height of stockpiles.

Discussion on the Parameter Setting
The volume calculation relies on the triangular mesh, the setting of the grid size for voxelization will have a certain influence on the experimental results. Taking real data as an example, grid size of voxelization was set with 0.6 m, 0.7, 0.8 m, 0.9 m, 1.0 m, 1.1 m, and 1.2 m respectively. The average PoV of different stockpiles with respect to different grid sizes are calculated. Because the number of point clouds obtained with different grid size are different, the IoU of the extraction results are no longer reported.
The calculation results are shown in Figure 28. With the decrease of grid size, the measurement error PoV in 28a decreases, but the time consumed by the algorithm increases in Figure 28b, It can be seen from Figure 28a that the PoV error does not change significantly when the grid size is around 0.7 m. But the time consumption will be greatly increased if the grid size becomes smaller than 0.7 m. Therefore, we suggest setting 0.7 m as the grid size for practical use. Actually, there are around 20 points per voxel if grid size equals to 0.7 m, which ensures the necessary point density for the proposed algorithm.

Conclusions
With the rapid increase of power supply demand, a large amount of stockpiles of coal have been formed during the process of coal excavation and transportation between the mines, ports, power plants and etc. Quantitative parameters, especially the volume of stockpile are important for the planning of coal production and consumption. Although laser scanning can collect the dense 3D coordinates of the stockpile surface for its quantification, stockpiles of coal have irregular shape, size, height, and base conditions, and may be overlapped with each other, which makes it hard to measure the different stockpiles automatically and accurately.
This paper proposes an algorithm to extract and measure the stockpiles from the 3D point cloud data using the multi-scale directional curvature. Although a stockpile has three parts of crest, slope, and base, it is distinctive and perceptible for the crest, which can then act as the starting point to extract and measure the stockpile. Firstly, the second-order directional curvature along multiple directions and at multiple scales is analyzed to extract the distinctive ridge of crest in the point cloud of stockpiles, resulting in invariant to the height, shape, size, and base conditions of different stockpiles. Then, starting with the crest points, a competitive growing strategy is proposed to accurately locate the points of slope in the stockpile, no matter whether there exist flatland, overlapping on the slope or not. Finally, the stockpile's volume is calculated by reconstructing the crest and slope points with many meshes and triangular prisms through the subsurface fitting and surface reconstruction, which is simple yet general enough for different stockpiles.
Through building a simulator system using a sandbox, Kinect depth sensor and projector, the effectiveness and accuracy of the proposed algorithm is comprehensively verified. Experiments on both the synthetic and real point cloud of stockpiles demonstrate that the proposed algorithm can extract the stockpiles with the average accuracy over 93.5% and measure the volume of stockpiles with the average precision over 93.7%. It is promising for automatically measuring the stockpiles like sand, soybean, etc., and facilitating the scientific storage and transportation management, as well as the maintenance of safety inventory during operation.
Although the proposed algorithm demonstrates its potential for extracting and measuring the stockpiles from point cloud, we suggest incorporating the analysis of boundary shape into the competitive strategy so that the boundary of a stockpile is not distorted a lot due to the over-growing.