A bulk material volume measurement method based on image features and 3D point cloud technology

By using image features and three-dimensional point cloud technology at the project site, combined with lidar and digital cameras, the rapid and accurate measurement of bulk material volume is achieved, and the problem of difficult to accurately count bulk material volume on the project site is solved, and the project quality and efficiency are improved .

CN114419130BActive Publication Date: 2024-12-06中国水利水电第七工程局有限公司
View PDF 0 Cites 2 Cited by

Patent Information

Application Number
CN202111578028.X
Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
Filing Date
2021-12-22
Publication Date
2024-12-06
Estimated Expiration
2041-12-22

AI Technical Summary

Technical Problem

At the project site, especially in the deep cover layer's ultra-standard vibration pile formation process, it is difficult to accurately count the bulk volume, resulting in difficulty in quality control and insufficient material use.

Method used

Using bulk volume measurement method based on image features and three-dimensional point cloud technology, the reference plane is obtained and triangulated by fixing three-dimensional lidar and digital camera on the forklift. The surface point cloud of bulk material is calculated and the volume of inclined triangular prism is realized to achieve fast and accurate measurement of bulk material volume.

Benefits of technology

It realizes rapid and accurate measurement of bulk material volume in complex engineering environments, improves project implementation efficiency and quality, and reduces the problem of insufficient material caused by human negligence.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN114419130B_ABST
    Figure CN114419130B_ABST
Patent Text Reader

Abstract

The present invention discloses a bulk material volume measurement method based on image features and three-dimensional point clouds, including: fixing a three-dimensional laser radar and a digital camera on a forklift to obtain the overall information of the bucket and its load; collecting and preprocessing point cloud data, optimizing the collected point cloud data through filtering and clustering algorithms; fusing image and point cloud information, calibrating parameters between the laser radar and the digital camera, projecting the point cloud to a pixel plane, fusing the straight line features of the image and the point cloud information; obtaining a reference plane; triangulating the bulk material surface point cloud, dividing the downsampled point cloud into triangular grids through greedy triangulation; calculating the volume: calculating the volume of the oblique truncated triangular prism corresponding to each triangular grid, and then accumulating to obtain the final bulk material volume. The method of the present invention has fast calculation speed, simple steps and accurate results; by collecting the point cloud and image information of the bucket, real-time volume calculation is realized, which is suitable for real-time monitoring in complex industrial and mining environments such as construction sites.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] The invention belongs to the technical field of engineering measurement, and in particular relates to a bulk material volume measurement method based on image features and three-dimensional point cloud technology. Background Art

[0002] At present, it is difficult to control the quality of the super-standard vibratory pile-forming process in the deep overburden on the project site, especially the deep pile-forming (35M), strong vibration (vibrator over 200KW), and irregular hole formation, which makes it difficult to accurately count the random filling volume on site. Therefore, a fast and accurate bulk material volume measurement method is needed. By rotating and extracting the bulk material volume, on the one hand, it can improve the implementation efficiency of the project, and on the other hand, it can avoid the problem of insufficient materials caused by human negligence, and improve the quality and speed of project construction.

[0003] In recent years, with the rapid development of 3D LiDAR technology, this sensor that can quickly obtain information about the surrounding environment has been widely used in many fields. As a sensor that actively obtains environmental information, LiDAR has better anti-interference ability and higher accuracy than traditional visual sensors, and can obtain depth information of the environment more quickly and accurately. Summary of the invention

[0004] The problem to be solved by the present invention is to meet the demand for fast and accurate measurement of bulk material volume during actual engineering implementation, and propose a volume measurement method based on image features and three-dimensional point cloud technology, which has fast calculation speed, accurate results, and can adapt to complex working environments such as construction sites.

[0005] The bulk material volume measurement method of the present invention is applied to bulk material backfilling or short-distance bulk material transfer construction in construction projects. Usually, bulk material transportation is carried out by a bucket excavator, and its shoveling and dumping are both mechanical operations. The present invention obtains fast and accurate measurement of bulk material volume through real-time measurement of the excavator bucket.

[0006] The present invention is achieved through the following technical solutions:

[0007] A bulk material volume measurement method based on image features and three-dimensional point cloud includes the following methods:

[0008] Step 1: A three-dimensional laser radar and a digital camera are fixedly installed on the forklift to obtain full-view information of the bucket and its load;

[0009] Step 2: Collection and preprocessing of point cloud data: Optimize the collected point cloud data through filtering and clustering algorithms;

[0010] Step 3: Fusion of image and point cloud information: Calibrate the parameters between the laser radar and the digital camera, project the point cloud onto the pixel plane, and fuse the straight line features of the image with the point cloud information;

[0011] Step 4: Obtain the reference plane;

[0012] Step 5: Triangulate the bulk material surface point cloud: divide the downsampled point cloud into triangular meshes through greedy triangulation;

[0013] Step 6. Calculate the volume: Calculate the volume of the truncated triangular prism corresponding to each triangular mesh, and then add them up to get the final bulk volume.

[0014] Further, the collection and preprocessing of point cloud data in step 2 includes:

[0015] Step 201: Considering that this method is continuously performed during the forklift operation, the acquisition time of the point cloud data is shortened to reduce the error caused by shaking. The acquisition time is set to 0.2 seconds, that is, the point cloud collected within 0.2 seconds is stored as one frame;

[0016] Step 202, the left and right width of the bucket is W. According to the movement characteristics of the forklift bucket, which only moves up and down and flips along the left and right axes, but does not translate left and right, the highest height that the bucket can reach is H, with the installation position of the radar as the origin. max , the lowest height that can be reached is H min , the closest distance to the radar is D max , the closest distance to the radar is D min , and use this to get the parameter X of the direct-pass filter min =D Min , X max =D max , Y min =-W / 2,Y max =W / 2, Z min =H min , Z max =H max ;The point cloud after straight-through filtering includes the bucket point cloud and part of the ground point cloud;

[0017] Step 203: using the Euclidean distance clustering algorithm, the point cloud is divided into two parts: the bucket and the ground, and the ground point cloud is removed, leaving the bucket bulk material point cloud data as the optimized collected point cloud data.

[0018] Further, the step 3 of fusing the image and point cloud information includes:

[0019] Step 301: according to the shape of the bucket, extract two adjacent boundaries of the bucket to determine the reference plane position; the two boundaries have a common intersection point O;

[0020] Calibrate the relationship between the laser radar and digital camera coordinate systems, and fix the relative positions of the laser radar and digital camera; the conversion relationship between the pixel coordinate system and the laser radar coordinate system is as follows:

[0021]

[0022] Where: f u and f v are the effective focal lengths in the horizontal and vertical directions, u0 and v0 are the coordinates of the center point of the image, R and T are the rotation matrix and translation matrix between the digital camera and the laser radar, respectively, which constitute the external parameter matrix of the camera; u and v are the pixel coordinates of the image; m ij The internal and external parameters are jointly transformed into matrix parameters to form matrix M;

[0023] Collect the point cloud and image information of the bucket. According to the relative projection relationship between the point cloud and image data, we can know that:

[0024]

[0025]

[0026] Where u and v are the pixel coordinates of the image, M is the transformation matrix obtained by calibration, and A is the target point in the point cloud that needs to be projected to the pixel plane.

[0027] Further, the step 4 of obtaining the reference plane includes:

[0028] Step 401, extract the bucket boundary in the image information through the LSD algorithm, and connect the line segments into a line; find the points that meet the projection to the straight line through the relative transformation relationship between the point cloud and the image; according to the imaging principle, the points projected from the point cloud to the same pixel on the image are roughly distributed in a line in space; fit the straight line through RANSAC and least squares, and keep the intersection of the two extracted straight lines;

[0029] Step 402: The collected bucket bulk material point cloud data does not exceed the boundary of the bucket, and the boundary of the point cloud is obtained by longitude and latitude scanning method. The specific process is as follows:

[0030] Project the point cloud onto the xoy plane, from X min To X max , with a step size of dx, find [x n ,x n +dx] in the interval min and max , save the index of the found point into the linked list; then min To Y max , with a step size of dy, find the x in the interval [y1,y1+dy] in turn min and x max , save the index of the found point into the linked list, and remove the index of the duplicate point;

[0031] Then project the point cloud data onto the yoz, xoz plane, repeat the above steps to find all the boundary points; perform multi-line fitting on the boundary point cloud, and take the two lines with the most line points as the reference lines; at the same time, determine the normal vector and parametric equation of the reference plane as Ax+By+Cz+D=0;

[0032] Step 403 : verify the obtained reference straight line through the conversion relationship of the camera radar coordinates in step 301 .

[0033] Further, the step 5 of triangulating the bulk material surface point cloud includes:

[0034] Step 501, dot product of the points in the point cloud and the normalized reference line, and the points whose dot product result is larger than the length and width of the bucket are out of range points and need to be removed;

[0035] Greedy triangulate the remaining points to obtain the information of which triangle the triangulated points belong to;

[0036] For the broken surface, voxel filtering and interpolation are performed respectively, where the voxel grid is set to 0.02*0.02*0.02 to obtain a surface with fewer holes;

[0037] Step 502, interpolation algorithm for broken faces: the three edges that make up a broken face belong to only one triangular face, while the three edges that make up a complete face belong to two triangular faces; find the points that make up each broken face; and then perform greedy triangulation on the points that make up the broken face.

[0038] Further, the step 6 of calculating the volume includes:

[0039] Step 601, calculate the volume of each oblique truncated triangular prism; take a triangular face ABC, the coordinates of A, B, and C are (X A ,Y A ,Z A )、(X B ,Y B ,Z B )、(X C ,Y C ,Z C ); The volume of the truncated triangular prism is obtained by adding the triangular prism AED-A'B'C' and the quadrangular pyramid A-BCED;

[0040] Add up the volumes of all the obliquely truncated triangular prisms and then add the volume of the standard bucket model to get the final bulk material volume.

[0041] In step 201 of the present invention, a Livox-Horizon laser radar is used to non-repetitively scan and collect a range of 81.7°×25.1° in front; the point cloud coordinate system takes the forward direction as the positive direction of the X-axis, the left direction as the positive direction of the Y-axis, and the upward direction as the positive direction of the Z-axis; the radar is placed so that the Y-axis of the radar coordinate system and the forklift are horizontally overlapped.

[0042] In step 203 of the present invention, for radial stray points generated at the boundary during the laser radar scanning process, the stray points are first removed by Gaussian statistical filtering, and then Euclidean clustering segmentation is performed based on connectivity in distance.

[0043] The specific steps of calibrating the parameters between the laser radar and the digital camera in step 301 of the present invention are: fixing the relative positions of the laser radar and the digital camera, collecting multiple sets of laser radar and digital camera data using the calibration board, manually selecting the corresponding point pair coordinates in the laser radar and the digital camera, and obtaining the conversion relationship between the radar point coordinates (x, y, z) and the camera pixel coordinates (u, v) through hyperstatic settlement, which is represented by the matrix M. The manually selecting the corresponding point pair coordinates in the laser radar and the digital camera is to select the corner point coordinates.

[0044] The bulk material volume measurement method based on image features of the present invention has fast calculation speed, simple steps and accurate results. By collecting point cloud and image information of the bucket, real-time volume calculation is realized, and the hardware calculation pressure is small. It is suitable for real-time volume measurement and progress monitoring in complex industrial and mining environments such as construction sites. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 Flow chart of the measuring method of the present invention;

[0046] Figure 2 The measurement method of the present invention extracts the schematic diagram of the edge of the bucket point cloud;

[0047] Figure 3 Schematic diagram of triangulation of surface point cloud of the measurement method of the present invention;

[0048] Figure 4 Schematic diagram of a truncated triangular prism according to the measurement method of the present invention. DETAILED DESCRIPTION

[0049] The present invention is further described below in conjunction with specific implementation methods. The specific implementation methods are further descriptions of the principles of the present invention and do not limit the present invention in any way. Technologies that are the same or similar to the present invention do not exceed the scope of protection of the present invention.

[0050] Combined with the attached pictures.

[0051] A bulk material volume measurement method based on image features and three-dimensional point cloud technology comprises the following steps:

[0052] Step 1: A three-dimensional laser radar and a digital camera are fixed on the forklift to obtain full-view information of the bucket and its load; the three-dimensional laser radar and the digital camera are fixed above the forklift, and the fixed position is preferably such that the full view of the bucket can be scanned when the bucket is naturally lifted.

[0053] Step 2: Collection and preprocessing of point cloud data; combined with Figure 2 ;

[0054] Step 201: Considering that this method is continuously carried out during the forklift operation, the acquisition time of the point cloud data is shortened to reduce the error caused by shaking. The acquisition time is set to 0.2 seconds, that is, the point cloud collected within 0.2 seconds is stored as one frame;

[0055] Step 202, the left and right width of the bucket is W. According to the movement characteristics of the forklift bucket, which only moves up and down and flips along the left and right axes, but does not translate left and right, the highest height that the bucket can reach is H, with the installation position of the radar as the origin. max , the lowest height that can be reached is H min , the closest distance to the radar is D max , the closest distance to the radar is D min , and use this to get the parameter X of the direct-pass filter min =D Min , X max =D max , Y min =-W / 2,Y max =W / 2, Z min =H min , Z max =H max The point cloud after straight-through filtering includes the bucket point cloud and a part of the ground point cloud.

[0056] Step 203: Since the bucket and bulk material point clouds are not connected to the ground point cloud, the point cloud can be divided into two parts, the bucket and the ground, by using the Euclidean distance clustering algorithm, and the ground point cloud is removed, leaving the bucket bulk material point cloud data for further processing.

[0057] Step 3: Fusion of image and point cloud information;

[0058] Step 301: Based on the shape of the bucket, extract two adjacent boundaries of the bucket to determine the reference plane position. The two boundaries have a common intersection point O.

[0059] Calibrate the relationship between the laser radar and digital camera coordinate systems, and fix the relative positions of the laser radar and digital camera. The conversion relationship between the pixel coordinate system and the laser radar coordinate system is as follows:

[0060]

[0061] Where: fu and f v are the effective focal lengths in the horizontal and vertical directions, u0 and v0 are the coordinates of the center point of the image, R and T are the rotation matrix and translation matrix between the digital camera and the laser radar, respectively, which constitute the external parameter matrix of the digital camera; u and v are the pixel coordinates of the image; m ij The internal and external parameters are jointly transformed into the matrix parameters which together form the matrix M.

[0062] Collect the point cloud and image information of the bucket. According to the relative projection relationship between the point cloud and image data, we can know that:

[0063]

[0064]

[0065] Where u and v are the pixel coordinates of the image, M is the transformation matrix obtained by calibration, and A is the target point in the point cloud that needs to be projected to the pixel plane.

[0066] Step 4: Obtain the reference plane;

[0067] Step 401, extract the bucket boundary in the image through the LSD algorithm, and connect the line segments into a line. Through the relative transformation relationship between the point cloud and the image, find the points that meet the projection to the straight line. According to the imaging principle, the points projected from the point cloud to the same pixel on the image are roughly distributed in a line in space. And use RANSAC and least squares to fit the straight line, and keep the intersection of the two extracted straight lines.

[0068] Step 402: The collected bucket bulk material point cloud data will not exceed the boundary of the bucket, so the boundary of the point cloud is obtained by longitude and latitude scanning method. The specific process is as follows:

[0069] Project the point cloud onto the xoy plane, from X min To X max , with a step size of dx, find [x n ,x n +dx] in the interval min and max , save the index of the found point into the linked list; then min To Y max , with a step size of dy, find the x in the interval [y1,y1+dy] in turn min and x max , save the index of the found point into the linked list, and remove the index of the duplicate point.

[0070] Then project the point cloud data onto the yoz, xoz planes, repeat the above steps to find all the boundary points. Perform multi-line fitting on the boundary point cloud, and take the two lines with the most line points as the reference lines. At the same time, the normal vector and parametric equation of the reference plane are determined, Ax+By+Cz+D=0.

[0071] Step 403 : verify the obtained reference straight line through the conversion relationship of the camera radar coordinates in step 301 .

[0072] Step 5: Triangulate the surface point cloud of bulk materials; Combine Figure 3 ;

[0073] Step 501: perform dot product of the points in the point cloud and the normalized reference straight line. Points whose dot product results are larger than the length and width of the bucket are out-of-range points and need to be removed.

[0074] Greedy triangulation is performed on the remaining points to obtain the information of which triangle the triangulated points belong to.

[0075] For the case of broken surfaces, voxel filtering and interpolation are performed respectively. Setting the voxel grid to 0.02*0.02*0.02 can obtain a surface with fewer holes.

[0076] Step 502, interpolation algorithm of broken surface:

[0077] The three edges that make up the broken face belong to only one triangle, while the three edges that make up the complete face belong to two triangles. In this way, the points that make up each broken face can be found. Then greedy triangle division is performed on the points that make up the broken face.

[0078] Step 6: Calculate the volume;

[0079] Step 601: Calculate the volume of each obliquely truncated triangular prism. Figure 4 As shown, take a triangular face ABC, the coordinates of A, B, and C are (X A ,Y A ,Z A )、(X B ,Y B ,Z B )、(X C ,Y C ,Z C ). The volume of the truncated triangular prism is obtained by adding the triangular prism AED-A'B'C' and the quadrangular pyramid A-BCED.

[0080] The two reference lines obtained in step 2 are OM and ON respectively, and the volume of the oblique triangular prism can be expressed as:

[0081]

[0082] in:

[0083]

[0084]

[0085]

[0086]

[0087] AF=2S A'B'C' / |B'C'|

[0088] Add up the volumes of all the truncated triangular prisms and then add the volume of the standard bucket model to get the final volume of the bulk material.

Claims

1. A bulk material volume measurement method based on image features and three-dimensional point clouds, characterized in that Includes the following methods: Step 1: A three-dimensional laser radar and a digital camera are fixedly installed on the forklift to obtain full-view information of the bucket and its load; Step 2: Collection and preprocessing of point cloud data: Optimizing the collected point cloud data through filtering and clustering algorithms; specifically including step 201, step 202, and step 203; Step 201, the duration of collecting information by the three-dimensional laser radar is 0.2 seconds; Step 202, the left and right width of the bucket is W. According to the movement characteristics of the forklift bucket, which only moves up and down and flips along the left and right axes but does not translate left and right, the highest height that the bucket can reach is H, with the installation position of the radar as the origin. max , the lowest height that can be reached is H min , the closest distance to the radar is D max , the closest distance to the radar is D min , and use this to get the parameter X of the direct-pass filter min =D Min , X max =D max , Y min =-W / 2,Y max =W / 2, Z min =H min , Z max =H max ; The point cloud after straight-through filtering includes the bucket point cloud and part of the ground point cloud; Step 203: using the Euclidean distance clustering algorithm, the point cloud is divided into two parts: the bucket and the ground, and the ground point cloud is removed, leaving the bucket bulk material point cloud data as the optimized collected point cloud data; Step 3: Fusion of image and point cloud information: calibrating the parameters between the laser radar and the digital camera, projecting the point cloud onto the pixel plane, and fusing the straight line features of the image with the point cloud information; specifically, step 301: Step 301: extract two adjacent boundaries of the bucket to determine the reference plane position according to the shape of the bucket; the two boundaries have a common intersection point O; Calibrate the relationship between the laser radar and digital camera coordinate systems, and fix the relative positions of the laser radar and digital camera; the conversion relationship between the pixel coordinate system and the laser radar coordinate system is as follows: Where: f u and f v are the effective focal lengths in the horizontal and vertical directions, u0 and v0 are the coordinates of the center point of the image, R and T are the rotation matrix and translation matrix between the digital camera and the laser radar, respectively, which constitute the external parameter matrix of the camera; u and v are the pixel coordinates of the image; m ij The internal and external parameters are jointly transformed into matrix parameters to form matrix M; Collect the point cloud and image information of the bucket. According to the relative projection relationship between the point cloud and image data, we can know that: Where u and v are the pixel coordinates of the image, M is the transformation matrix obtained by calibration, and A is the target point in the point cloud that needs to be projected to the pixel plane; Step 4: Obtaining a reference plane; specifically including step 401, step 402, and step 403: Step 401, extract the bucket boundary in the image information through the LSD algorithm, and connect the line segments into a line; find the points that meet the projection to the straight line through the relative transformation relationship between the point cloud and the image; according to the imaging principle, the points projected from the point cloud to the same pixel on the image are roughly distributed in a line in space; fit the straight line through RANSAC and least squares, and keep the intersection of the two extracted straight lines; Step 402: The collected bucket bulk material point cloud data does not exceed the boundary of the bucket, and the boundary of the point cloud is obtained by longitude and latitude scanning method. The specific process is as follows: Project the point cloud onto the xoy plane, from X min To X max , with a step size of dx, find [x n ,x n +dx] in the interval min and max , save the index of the found point into the linked list; then min To Y max , with a step size of dy, find the x in the interval [y1,y1+dy] in turn min and x max , save the index of the found point into the linked list, and remove the index of the duplicate point; Then project the point cloud data onto the yoz, xoz plane, repeat the above steps to find all the boundary points; perform multi-line fitting on the boundary point cloud, and take the two lines with the most line points as the reference lines; at the same time, determine the normal vector and parametric equation of the reference plane as Ax+By+Cz+D=0; Step 403: verify the obtained reference straight line through the conversion relationship of the camera radar coordinates in step 301; Step 5: Triangulate the bulk material surface point cloud: divide the downsampled point cloud into triangular meshes through greedy triangulation; Step 6. Calculate the volume: Calculate the volume of the truncated triangular prism corresponding to each triangular mesh, and then add them up to get the final bulk volume.

2. The bulk material volume measurement method based on image features and three-dimensional point cloud according to claim 1 is characterized in that The triangulated bulk material surface point cloud in step 5 includes: Step 501, dot product of the points in the point cloud and the normalized reference line, and the points whose dot product result is larger than the length and width of the bucket are out of range points and need to be removed; Greedy triangulate the remaining points to obtain the information of which triangle the triangulated points belong to; For the broken surface, voxel filtering and interpolation are performed respectively, where the voxel grid is set to 0.02*0.02*0.02 to obtain a surface with fewer holes; Step 502, interpolation algorithm for broken faces: the three edges that make up a broken face belong to only one triangular face, while the three edges that make up a complete face belong to two triangular faces; find the points that make up each broken face; and then perform greedy triangulation on the points that make up the broken face.

3. The bulk material volume measurement method based on image features and three-dimensional point cloud according to claim 2 is characterized in that The volume calculation described in step 6 includes: Step 601, calculate the volume of each oblique truncated triangular prism; take a triangular face ABC, the coordinates of A, B, and C are (X A ,Y A ,Z A )、(X B ,Y B ,Z B )、(X C ,Y C ,Z C ); The volume of the truncated triangular prism is obtained by adding the triangular prism AED-A'B'C' and the quadrangular pyramid A-BCED; Add up the volumes of all the truncated triangular prisms and then add the volume of the standard bucket model to get the final bulk material volume.

4. The bulk material volume measurement method based on image features and three-dimensional point clouds according to claim 3 is characterized in that: In step 201, a Livox-Horizon laser radar is used to non-repetitively scan and collect the range of 81.7°×25.1° in front; the point cloud coordinate system takes the forward direction as the positive direction of the X axis, the left direction as the positive direction of the Y axis, and the upward direction as the positive direction of the Z axis; the radar is placed so that the Y axis of the radar coordinate system and the forklift are horizontally overlapped.

5. The bulk material volume measurement method based on image features and three-dimensional point clouds according to claim 2, characterized in that: In step 203, for radial stray points generated at the boundary during the laser radar scanning process, the stray points are first removed by Gaussian statistical filtering, and then Euclidean clustering segmentation is performed based on connectivity in distance.

6. The bulk material volume measurement method based on image features and three-dimensional point clouds according to claim 2, characterized in that: The specific steps of calibrating the parameters between the laser radar and the digital camera in step 301 are as follows: The relative positions of the laser radar and digital camera are fixed, and multiple sets of laser radar and digital camera data are collected using the calibration plate. The corresponding point pair coordinates in the laser radar and digital camera are manually selected, and the conversion relationship between the radar point coordinates (x, y, z) and the camera pixel coordinates (u, v) is obtained through hyperstatic settlement, which is represented by the matrix M.

7. The bulk material volume measurement method based on image features and three-dimensional point clouds according to claim 6, characterized in that: The coordinates of the corner points are selected by manually selecting the corresponding point pairs in the laser radar and the digital camera.

Citation Information

Cited By

  • Shovel loading volume obtaining method and system based on material pile three-dimensional reconstruction before and after shovel loading

    CN116258832A

  • Method and System for Obtaining Loading Volume Based on 3D Reconstruction of Material Piles before and after Loading

    CN116258832B