An automatic measurement method for part phenotype parameters of pinus massoniana based on three-dimensional point cloud image
By using 3D point cloud image processing technology, the phenotypic parameters of Masson pine plants are automatically identified and calculated, solving the problems of time consumption and human error in traditional methods, and achieving efficient and accurate phenotypic detection of Masson pine.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- ZHEJIANG UNIV OF SCI & TECH
- Filing Date
- 2022-07-19
- Publication Date
- 2026-06-26
AI Technical Summary
In existing technologies, the three-dimensional phenotypic measurement methods for Masson pine are time-consuming, labor-intensive, and prone to human error. There is a lack of fast and accurate automated detection solutions, especially when needle-like leaves and tillers are obstructing the data, making it difficult to obtain high-quality and repeatable phenotypic data.
An automatic measurement method based on three-dimensional point cloud images is adopted. Through steps such as acquisition, preprocessing, registration, skeletonization and leaf separation, Masson pine plants are automatically identified and the main diameter, ground diameter, average leaf length and leaf number are calculated to achieve digital phenotypic evaluation.
It improves the efficiency and accuracy of phenotypic detection of Masson pine, provides rapid and low-cost non-destructive testing, is applicable to plants of different shapes, sizes and growth stages, and has better grading accuracy.
Smart Images
Figure CN115393418B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to a non-destructive and automated detection method for phenotypic parameters of Masson pine in the field of three-dimensional point cloud images, and in particular to an automatic measurement method for some phenotypic parameters of Masson pine based on three-dimensional point cloud images. Background Technology
[0002] Masson pine, also known as Chinese pine, mountain pine, and fir pine, is distributed in areas south of the Qinling Mountains and Huai River in my country. It is not only an important timber crop for industrial and agricultural production but also a plant with high medicinal value. Currently, selecting varieties and seedlings based on germplasm resources is a crucial means to improve the Masson pine industry. By analyzing the relationship between genotype and phenotype, high-yielding and stress-resistant Masson pine varieties can be rapidly cultivated. However, rapid phenotypic analysis methods for Masson pine (such as rapid and non-destructive measurements of plant height, diameter at breast height, and leaf area) remain a bottleneck in genetics. Traditional phenotypic analysis methods are often time-consuming, labor-intensive, and prone to human error. Therefore, developing fully automated technologies for collecting crop phenotypic data is essential, especially visual phenotyping technologies with higher accuracy and efficiency.
[0003] While two-dimensional (2D) images have been widely used in image-based phenotyping techniques, three-dimensional (3D) plant analysis offers greater stability and accuracy. Beyond qualitative description, 3D imaging can measure more features and traits, such as for 3D canopy reconstruction, LAI estimation, and the measurement of plant traits like height, volume, biomass, stem diameter, and leaf area. However, after reviewing domestic and international literature, no systematic 3D phenotyping scheme for Masson pine has been found. In 3D phenotyping of Masson pine, occlusion caused by needle-like leaves and tillers often necessitates more sophisticated methods. To obtain high-quality, repeatable phenotypic datasets, suitable collection areas need to be manually selected, followed by data standardization using analysis software to improve reliability. This process still requires personnel with certain biological and technical backgrounds. Therefore, manual measurements are generally limited in scale and throughput, and operators must follow standardized procedures to complete the collection task. Currently, there is a lack of a rapid phenotyping scheme for Masson pine to measure parameters such as main diameter. Summary of the Invention
[0004] To address the problems existing in the background technology, the purpose of this invention is to provide an automatic measurement method for some phenotypic parameters of Masson pine based on three-dimensional point cloud images. This method can automatically identify Masson pine plants in three-dimensional point cloud images and complete the calculation and classification of digital phenotypic evaluation indicators such as main diameter length, ground diameter, average leaf length, and leaf length. This improves detection efficiency and, together with imaging and other appearance detection methods, lays a technical foundation for the application of Masson pine phenomics.
[0005] The technical solution adopted in this invention includes the following steps:
[0006] 1) Collect three-dimensional point cloud images of Masson pine plants. During the collection process, rotate the Masson pine plants and record one three-dimensional point cloud image for each rotation. Perform N rotations on each Masson pine plant and record N three-dimensional point cloud images.
[0007] 2) Preprocess each 3D point cloud image to obtain the 3D point cloud of the Masson pine plant itself;
[0008] 3) Register the N preprocessed 3D point clouds to obtain a fused point cloud;
[0009] 4) Extract the minimum spanning tree of candidate points and the skeleton points of the main stem, and then obtain the point cloud of the main stem and the point cloud of the leaves;
[0010] 5) Obtain the main stem length H, diameter at ground level of Masson pine, average leaf length and number of leaves by processing the minimum spanning tree of candidate points, main stem point cloud and leaf point cloud.
[0011] In this method, the Masson pine plant is placed in the soil of a flowerpot, which is placed on a table with the soil just overflowing from the pot and the top surface of the soil level with the top surface of the flowerpot.
[0012] The preprocessing in step 2) specifically involves:
[0013] 2.1) Use a pass-through filter to remove the background from the 3D point cloud image;
[0014] 2.2) The RANSAC (Random Sample Consensus) plane fitting method is used to detect the desktop plane in the 3D point cloud image after removing the background, and then the desktop plane is removed by plane segmentation.
[0015] 2.3) A radius-based filtering method is used to remove a small number of sparsely distributed noise points in the 3D point cloud image after removing the background and the desktop plane;
[0016] 2.4) Obtain the lowest point of the pot where the Masson pine plant is located from the 3D point cloud image after removing noise points. Then, estimate the soil plane height of the pot based on the lowest point and the pot height to determine the soil plane. Take the area above the soil plane as the region of interest and remove the area below the soil plane using a pass-through filter. All points in the region of interest constitute the 3D point cloud of the Masson pine plant itself.
[0017] Step 3) specifically refers to:
[0018] 3.1) For both the source and target point clouds, traverse each point and perform the following processing to obtain the feature points of the source and target point clouds:
[0019] Points in the point cloud are selected as processing points. A region point cloud is formed by considering all points within a neighborhood radius of the processing point. Principal component analysis (PCA) is used to estimate the covariance matrix of the region point cloud. The eigenvector corresponding to the smallest eigenvalue of the covariance matrix is then used as the normal vector of the processing point. Next, the angles between the normal vector of the processing point and the normal vectors of all points within the region point cloud and its neighbors are calculated. The mean of all these angles is then calculated as the average angle. An angle threshold is set, and the following judgments are made:
[0020] If the mean angle between the processing points is greater than the angle threshold, the processing points are retained as feature points; the larger the angle between the normal vectors, the greater the curvature of the region; otherwise, the processing points are not retained as feature points.
[0021] Repeat step 3.1) to process all points in the point cloud and obtain all feature points.
[0022] 3.2) Calculate the feature descriptors of feature points in the source point cloud and the target point cloud, and use the feature descriptors in the source point cloud and the target point cloud to perform coarse registration using the sampling consistency method SAC-IA;
[0023] 3.3) The source point cloud and target point after coarse registration coordinate transformation are used as input, and the nearest point iterative algorithm ICP method is used for fine registration. In the process of iteratively searching for the nearest point in the ICP method, the KD-tree based neighbor search algorithm is used for improvement to obtain the transformation matrix between the source point cloud and the target point cloud.
[0024] 3.4) Take the last 3D point cloud in the N preprocessed 3D point clouds as the target point cloud, and the remaining 3D point clouds as the source point clouds. Traverse each source point cloud and repeat the above steps 3.1) to 3.4) to perform registration processing between the source point clouds and the target point clouds to obtain the transformation matrix between each pair of source point clouds and target point clouds.
[0025] 3.5) Register all source point clouds to the target point cloud using their respective transformation matrices to obtain each registered source point cloud. Combine the registered source point clouds and target point clouds into a complete fused point cloud by point cloud addition.
[0026] Step 4) specifically refers to:
[0027] 4.1) The fused point cloud is sliced along the growth direction for the first time. The number of slices in the growth direction is customized. Euclidean clustering is performed on the point cloud in each slice layer obtained in the first slice to obtain each cluster.
[0028] 4.2) Based on the clusters obtained in step 4.1), make the following judgments:
[0029] For clusters in step 4.1) where the number of points is less than or equal to the preset threshold, the centroid of the cluster is directly extracted as candidate points for the main stem skeleton.
[0030] For clusters with more than the preset threshold number of points obtained in step 4.1), a second slice is made along any direction perpendicular to the growth direction (X direction). The number of slices in the X direction is customized. The points in each slice layer obtained in the second slice are clustered according to the gray value to obtain each cluster. The average gray value of each cluster is then calculated. The centroid of the slice layer of the second slice containing the cluster with the largest gray value is taken as the candidate point of the main stem skeleton.
[0031] The set of candidate points is composed of all the candidate points of the main stem skeleton;
[0032] 4.3) Construct the minimum spanning tree of the candidate point set, and then use the DAG single-source longest path algorithm to search for the longest path in the minimum spanning tree of the candidate points. All candidate points of the longest path form a set of main stem skeleton points containing the canopy.
[0033] 4.4) Set a canopy length threshold. Search downwards from the top candidate point in the set of main stem skeleton points containing the canopy. When a candidate point is found, if the search path length is greater than the canopy length threshold, stop the search. All candidate points found in the longest path up to the top candidate point are discarded. The remaining candidate points in the longest path are used as the skeleton points of the main stem.
[0034] 4.5) Interpolate between the skeleton points of each main stem, inserting a point at fixed intervals; after completing the interpolation of the skeleton points of the main stem, use the points obtained by searching along the skeleton points from bottom to top using the radius search method of Kdtree as the points in the main stem point cloud, thereby obtaining the main stem point cloud;
[0035] The fixed distance is no greater than the estimated radius of the main stem at the bottom of the Masson pine at 0.02-0.05m.
[0036] 4.6) Subtract the main stem point cloud from the fused point cloud obtained after registration in step 3) to obtain the point cloud of all leaves, i.e., the leaf point cloud, thus realizing the separation of stem and leaf of Masson pine.
[0037] In step 5), the length H of the main stem is obtained as follows: Based on the skeleton points of the main stem, a minimum spanning tree is built from the skeleton points of the main stem as the minimum spanning tree of the main stem skeleton. The distance between each two adjacent skeleton points is extracted according to the point order in the minimum spanning tree of the main stem skeleton. The sum of the distances between all two adjacent skeleton points is obtained as the length H of the main stem.
[0038] In step 5), the diameter at ground level of the Masson pine is obtained in the following manner:
[0039] S11. Take the skeleton point at a preset fixed distance from the lowest skeleton point in the minimum spanning tree of the candidate points as the starting point, draw a line connecting the starting point and the skeleton point adjacent to the minimum spanning tree of the candidate points above it as the lower normal, and generate a plane perpendicular to the lower normal through the starting point as the lower plane.
[0040] S12. Take the nearest skeleton point along the minimum spanning tree of the candidate points above the starting point as the secondary starting point, and draw a line connecting the secondary starting point and the nearest skeleton point along the minimum spanning tree of the candidate points above it as the upper normal. Generate a plane perpendicular to the upper normal through the secondary starting point as the upper plane.
[0041] S13. In the main stem point cloud, all points in the main stem point cloud located between the upper and lower planes are vertically projected onto the lower plane to obtain projection points. The maximum distance between any two projection points is calculated as the ground diameter of the Masson pine.
[0042] In step 5), the average blade length is obtained in the following manner:
[0043] S21. Select a skeleton point from the minimum spanning tree of candidate points as the starting point, draw a line connecting the starting point and the skeleton point adjacent to the minimum spanning tree of candidate points above it as the lower normal, and generate a plane perpendicular to the lower normal through the starting point as the lower plane.
[0044] S22. Take the nearest skeleton point of the minimum spanning tree of the candidate points above the starting point as the secondary starting point, and draw a line connecting the secondary starting point and the nearest skeleton point of the minimum spanning tree of the candidate points above it as the upper normal. Generate a plane perpendicular to the upper normal through the secondary starting point as the upper plane.
[0045] S23. Project all points in the main stem point cloud located between the upper and lower planes vertically onto the lower plane to obtain projection points. Construct the smallest α-bounding box that encloses all projection points. Calculate the average distance between all points on the boundary of the smallest α-bounding box and the projection point (only one point) projected onto the lower plane from the starting point:
[0046] If the average distance is greater than the preset distance threshold, the average distance is used as the average leaf length of the region where the starting point is located.
[0047] If the average distance is less than or equal to the preset distance threshold, it is assumed that there are no blades in the area where the starting point is located.
[0048] S24. Repeat steps S11-S14 for each skeleton point in the minimum spanning tree of the candidate points to obtain the average leaf length of the region where each skeleton point is located.
[0049] In step 5), the number of blades is obtained in the following manner:
[0050] S31. Select a skeleton point from the minimum spanning tree of candidate points as the starting point, draw a line connecting the starting point and the skeleton point adjacent to the minimum spanning tree of candidate points above it as the lower normal, and generate a plane perpendicular to the lower normal through the starting point as the lower plane.
[0051] S32. Take the nearest skeleton point of the minimum spanning tree of the candidate points above the starting point as the secondary starting point, and draw a line connecting the secondary starting point and the nearest skeleton point of the minimum spanning tree of the candidate points above it as the upper normal. Generate a plane perpendicular to the upper normal through the secondary starting point as the upper plane.
[0052] S33. Count the number of points in the leaf point cloud located between the upper and lower planes.
[0053] If the number of points is greater than the preset number of points threshold, the average distance density of all points in the leaf point cloud located between the upper and lower planes is calculated, and the number of leaves in the starting area is obtained by multiplying the average distance density by a parameter.
[0054] If the number of points is not greater than the preset point threshold, then there are no leaves in the region where the starting point is located;
[0055] The average distance density is obtained by summing the minimum distances between each point and other points and dividing by the number of points.
[0056] S34. Repeat steps S11-S14 for each skeleton point in the minimum spanning tree of the candidate points to obtain the number of leaves in the region of each skeleton point.
[0057] The number of blades in the region of each skeleton point is calculated using the above method, and the results are summed to obtain the total number of blades.
[0058] The beneficial effects of this invention are:
[0059] Masson pine is one of my country's important afforestation plants with high economic value. The method of this invention helps improve the efficiency of cultivating superior Masson pine seedlings. By analyzing the growth of plants grown in different environments, high-yielding and stress-resistant plants can be cultivated more quickly. Currently, plant growth is mainly reflected by phenotypic parameters such as plant height and ground diameter, which are originally obtained primarily through phenotypic measurement methods. Traditional phenotypic measurement operations are time-consuming, labor-intensive, and prone to human error.
[0060] This invention utilizes a camera to obtain three-dimensional point cloud images of plants, enabling high-precision and high-efficiency crop phenotypic data acquisition, and achieving rapid, low-cost, accurate, and non-destructive detection of crop phenotypic parameters.
[0061] The method of this invention adopts an automated point cloud image processing scheme, including preprocessing, registration, skeletonization, main stem extraction, leaf extraction and other processes. After that, multiple phenotypic parameters are automatically calculated and corresponding grading strategies are proposed. It has universality for Masson pine plants of different shapes, sizes and stages, and can provide digital evaluation indicators. It has better grading accuracy than other methods. Attached Figure Description
[0062] Figure 1 This is a flowchart of the method of the present invention.
[0063] Figure 2 In this embodiment of the invention, the image is a point cloud image of a Masson pine tree.
[0064] Figure 3 This is an image showing the effect after the pretreatment step in an embodiment of the present invention. The soil part is marked with a darker color, indicating that it has actually been removed.
[0065] Figure 4 This is a frontal point cloud image of the registered plant in the embodiment;
[0066] Figure 5 This is a point cloud image of the main stem extracted in an embodiment of the present invention;
[0067] Figure 6 This is a point cloud image of the blade extracted in an embodiment of the present invention. Detailed Implementation
[0068] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
[0069] like Figure 1 As shown, embodiments of the present invention are as follows:
[0070] 1) Collect three-dimensional point cloud images of Masson pine plants. During the collection process, rotate the Masson pine plants and record one three-dimensional point cloud image for each rotation. Perform N rotations on each Masson pine plant and record N three-dimensional point cloud images.
[0071] In this embodiment, the Masson pine plant is placed in the soil of a flowerpot, which is placed on a table with the soil just overflowing from the pot and the top surface of the soil level with the top surface of the flowerpot.
[0072] A 3D plant model was created using a Kinect camera and a manual turntable under the following conditions: The camera position was fixed, and the plant was placed on the turntable. The Kinect height was set to 0.7 meters. The camera angle on the Yunteng 691 support was adjusted to a downward-looking angle of 20 degrees (i.e., a vertical angle of 70 degrees with the world coordinate system). The distance between the Kinect and the plant was approximately 0.5-0.6 meters, ensuring the entire plant was centered in the camera's field of view, with appropriate adjustments made for different plants. The distance between the turntable and the screen was 1.8-3 meters, ensuring the plant's silhouette was entirely within the screen's area, minimizing background interference. The plant was placed on the turntable, and a point cloud image was captured every 60 degrees of rotation, resulting in 6 point cloud blocks. The number of point clouds could be adjusted as needed, and the point cloud was captured using the SDK included with the DK camera.
[0073] 2) Preprocess each 3D point cloud image to obtain the 3D point cloud of the Masson pine plant itself;
[0074] 2.1) Use a pass-through filter to remove the background from the 3D point cloud image;
[0075] The main object is the plant, but the scene includes a backdrop that serves no purpose for experimental measurements and interferes with subsequent point cloud analysis of the plant. A direct-pass filter is used to remove the backdrop.
[0076] The input parameters for the pass-through filter are a 3D point cloud image. The pass-through parameters are as follows: d is the distance between the camera and the screen, and x is the distance between the camera and the plants. The filtering range is (dx) / d. In this experiment, d is 1.8 and x is 0.5.
[0077] 2.2) The RANSAC (Random Sample Consensus) plane fitting method is used to detect the desktop plane in the 3D point cloud image after removing the background, and then the desktop plane is removed by plane segmentation.
[0078] The experimental platform and table are measurement-irrelevant variables, having no effect on the measurement process and interfering with subsequent point cloud analysis of the plants; therefore, they need to be removed. The experimental platform (i.e., the table) was located in the plane point cloud using RANSAC (a plane fitting method), and thus removed. The target geometry was set to SACMODEL_PLANE, and the segmentation method was set to random sampling. The maximum number of iterations and the threshold range were set to their default values, and the resulting plane was selected to segment the remaining points.
[0079] 2.3) A radius-based filtering method is used to remove a small number of sparsely distributed noise points in the 3D point cloud image after removing the background and the desktop plane;
[0080] After removing the background and experimental platform, there will be some discrete noise. These discrete point clouds are invalid data and will interfere with Euclidean clustering in the subsequent skeletonization. Here, we use radius filtering to remove these discrete point clouds with a radius range of 0.5 and a search neighborhood of 100.
[0081] 2.4) Obtain the lowest point of the pot where the Masson pine plant is located from the 3D point cloud image after removing noise points. Then, estimate the soil plane height of the pot based on the lowest point and the pot height to determine the soil plane. Take the area above the soil plane as the region of interest and remove the area below the soil plane using a pass-through filter. All points in the region of interest constitute the 3D point cloud of the Masson pine plant itself.
[0082] 3) Register the N preprocessed point cloud images;
[0083] 3.1) For both the source and target point clouds, traverse each point and perform the following processing to obtain the feature points of the source and target point clouds:
[0084] Points in the point cloud are selected as processing points. A region point cloud is formed by considering all points within a neighborhood radius of the processing point. Principal component analysis (PCA) is used to estimate the covariance matrix of the region point cloud. The eigenvector corresponding to the smallest eigenvalue of the covariance matrix is then used as the normal vector of the processing point. Next, the angles between the normal vector of the processing point and the normal vectors of all points within the region point cloud and its neighbors are calculated. The mean of all these angles is then calculated as the average angle. An angle threshold is set, and the following judgments are made:
[0085] If the mean angle between the processing points is greater than the angle threshold, the processing points are retained as feature points; the larger the angle between the normal vectors, the greater the curvature of the region; otherwise, the processing points are not retained as feature points.
[0086] Repeat step 3.1) to process all points in the point cloud and obtain all feature points.
[0087] Specifically, before the initial registration, due to the large amount of point cloud data, the registration efficiency is reduced. The original point cloud is downsampled to effectively reduce the number of points and improve the registration efficiency.
[0088] Calculate the normal vector of the downsampled point cloud: Set up a point cloud normal vector calculation unit, set the search method of the calculation unit to kdTree, set the search range and the number of points to be searched. In this experiment, the number of points is set to 10. Calculate the normal vector of the 10 points searched by kdTree within the search range.
[0089] Calculate the angle and mean of the normal vectors, and determine whether they are feature points. If the angle of the point cloud normal vectors in a local area is large, it indicates that the area is relatively undulating; conversely, if the angle of the normal vectors does not change much, it indicates that the area is relatively flat.
[0090] 3.2) Calculate the FPFH feature descriptors for feature points in the source and target point clouds;
[0091] After extracting feature points based on the angle characteristics of the point cloud normal vectors, FPFH is used to describe these feature points. FPFH is a method for describing local features by extracting the neighborhood geometry of feature points. It constructs a feature histogram using the normal vectors of the point cloud and the surface curvature of neighboring points, thereby quickly determining the correspondence between feature points.
[0092] This uses OMP acceleration. OpenMP's Parallelfor instruction is a standard fork / join parallel mode. The basic idea is that the program starts with only one main thread, and the serial parts of the program are executed by the main thread. Parallel parts are executed by spawning other threads, but the serial parts will not be executed until the parallel parts have finished. In other words, the non-parallel parts can only be executed after all the parallel parts of the OpenMP program have finished. This speeds up the feature point extraction process in FPFH.
[0093] 3.3) Using the feature descriptors in the source and target point clouds, coarse registration is performed using the Sampling Consistency Method (SAC-IA).
[0094] (1) Selecting sampling points
[0095] Select n sampling points from the source point cloud P. In order to ensure that the sampling points have different FPFH features, the distance between the sampling points should be greater than the minimum distance threshold d.
[0096] (2) Corresponding point search
[0097] Find one or more points in the target point cloud Q that have similar PFFH features to the point sampling point, and use these similar points as the corresponding points of the sampling point in the point cloud Q.
[0098] (3) Calculate the rigid body transformation matrix
[0099] The transformation matrix is calculated based on the corresponding points. The registration performance is judged by the transformation parameter "sum of distance errors" function, and the Huber formula is used to determine the "sum of distance errors" function. The optimal solution is selected based on the transformation of the corresponding points to minimize the point cloud registration distance error, and the final rigid body transformation matrix is obtained, completing the coarse registration.
[0100] 3.4) The source point cloud and target point after coarse registration coordinate transformation are used as input, and fine registration is performed using the ICP method. In the iterative search for the nearest point in the ICP method, the KD-tree based proximity search algorithm is used for improvement to obtain the transformation matrix between the source point cloud and the target point cloud.
[0101] This invention optimizes the traditional ICP algorithm for accurate registration by using a K-dimensional tree nearest neighbor search algorithm, which accelerates the search for corresponding point pairs and effectively improves computational efficiency.
[0102] 3.5) Take the last 3D point cloud in the N preprocessed 3D point clouds as the target point cloud, and the remaining 3D point clouds as the source point clouds. Traverse each source point cloud and repeat the above steps 3.1) to 3.4) to perform registration processing between the source point clouds and the target point clouds to obtain the transformation matrix between each pair of source point clouds and target point clouds.
[0103] 3.6) Register all source point clouds to the target point cloud using their respective transformation matrices to obtain each registered source point cloud. Combine the registered source point clouds and target point clouds into a complete fused point cloud by point cloud addition.
[0104] 4) Based on the idea of skeletonization, extract the skeleton points of the main stem and the skeleton points of the leaves; then obtain the point cloud of the main stem and the leaves;
[0105] 4.1) The fused point cloud is sliced along the growth direction for the first time. The number of slices in the growth direction is customized. Euclidean clustering is performed on the point cloud in each slice layer obtained in the first slice to obtain each cluster.
[0106] 4.2) Based on the clusters obtained in step 4.1), make the following judgments:
[0107] For clusters in step 4.1) where the number of points is less than or equal to the preset threshold, the centroid of the cluster is directly extracted as candidate points for the main stem skeleton.
[0108] For clusters with more than the preset threshold number of points obtained in step 4.1), a second slice is made along any direction perpendicular to the growth direction (X direction). The number of slices in the X direction is customized. The points in each slice layer obtained in the second slice are clustered according to the gray value to obtain each cluster. The average gray value of each cluster is then calculated. The centroid of the slice layer of the second slice containing the cluster with the largest gray value is taken as the candidate point of the main stem skeleton.
[0109] The set of candidate points is composed of all the candidate points of the main stem skeleton.
[0110] The main stem of the Masson pine has a unique characteristic: its point cloud intensity value is the highest. After slicing in the Y direction, each adjacent slice is then horizontally sliced in the X direction, and the intensity value of each horizontal slice is calculated. Here, the average value of the RGB three channels is used as the intensity value of each horizontal slice. The horizontal slice with the highest intensity value is the slice containing the main stem of the Masson pine. Finally, the centroid of this slice is calculated as the skeletal point of the main stem of the Masson pine.
[0111] 4.3) Construct the minimum spanning tree of the candidate point set, and then use the DAG single-source longest path algorithm to search for the longest path in the minimum spanning tree of the candidate points. All candidate points of the longest path form a set of main stem skeleton points containing the canopy.
[0112] The main stem of the Masson pine is characterized by its longest length. Based on this characteristic, we first construct an undirected single-source path graph starting from the bottom node. Since there are no closed loops in the graph, we use a Directed Acyclic Graph (DAG) for the critical path. We construct an array dp[n+1], where dp[i] represents the longest path starting from this node. To find dp[i], we only need to find dp[j] of the next node starting from i, such that dp[i] = max(dp[i], dp[j] + G[i][j]). Through this state transition equation, we solve for the longest single-source path. All the skeleton points on this longest path are the skeleton points of the main stem.
[0113] 4.4) Set a canopy length threshold. Search downwards from the top candidate point in the set of main stem skeleton points containing the canopy. When a candidate point is found, if the search path length is greater than the canopy length threshold, stop the search. In the longest path, all candidate points found up to the top candidate point are discarded. The remaining candidate points in the longest path are used as the skeleton points of the main stem.
[0114] The main stem of the Masson pine is characterized by its longest length. Based on this characteristic, we first construct an undirected single-source path graph starting from the bottom node. Since the graph contains no closed loops, we use a Directed Acyclic Graph (DAG) with critical paths. We create an array `dp[n+1]`, where `dp[i]` represents the longest path starting from this node. To find `dp[i]`, we only need to find `dp[j]` of the next node starting from `i`, such that `dp[i] = max(dp[i], dp[j] + G[i][j])`. Using this state transition equation, we find the longest single-source path. All the skeleton points on this longest path are the skeleton points of the main stem. We then traverse from the top of the main stem's skeleton points, deleting skeleton points within the length of the canopy at the top. The resulting skeleton points are the skeleton points of the main stem.
[0115] 4.5) Interpolate between the skeleton points of each main stem, inserting a point at fixed intervals, with the fixed interval not exceeding the estimated radius of the main stem at the bottom of the Masson pine at 0.02-0.05m; after completing the interpolation of the skeleton points of the main stem, use the points obtained by searching along the skeleton points from bottom to top using the radius search method of Kdtree as the points in the main stem point cloud, thereby obtaining the main stem point cloud;
[0116] After extracting the skeleton points of the main stem, the point cloud of the main stem of the Masson pine is restored based on the skeleton points of the main stem. The restoration method here is the radius search of Kdtree, which takes the radius of the main stem at the bottom of the Masson pine (0.02-0.05m) as the search range and searches for all points within the radius of the skeleton points.
[0117] However, this method has a drawback: the skeleton points of the main stem are relatively sparse, and the radius search using the kd-tree of the main stem skeleton points cannot find all the points of the main stem. Therefore, the skeleton points of the main stem need to be expanded. The expansion method is as follows: for the lines connecting each skeleton point, expand by one point at regular intervals, with the interval being the radius of the main stem at the base of the Masson pine (0.02-0.05m). After this expansion, the skeleton points of the main stem can be used to search the point cloud containing all the main stems.
[0118] 4.6) Subtract the main stem point cloud from the fused point cloud obtained after registration in step 3) to obtain the point cloud of all leaves, i.e., the leaf point cloud, thus realizing the separation of stem and leaf of Masson pine.
[0119] For the expanded main stem skeleton points, these skeleton points are used as the search center for radius search of kdtree. The main stem radius at the bottom of the Masson pine (0.02-0.05m) is used as the search range. The searched point cloud is extracted. These point clouds are the point clouds where the main stem is located. Then, the point cloud where the main stem is located is subtracted from the original point cloud to obtain the leaf point cloud, thus realizing the separation of the stem and leaf of Masson pine.
[0120] 5) Obtain the main stem length H, diameter at ground level of Masson pine, average leaf length and number of leaves by processing the minimum spanning tree of candidate points, main stem point cloud and leaf point cloud.
[0121] A. Main stem length H: Based on the skeleton points of the main stem, construct a minimum spanning tree from the skeleton points of the main stem as the minimum spanning tree of the main stem skeleton. Extract the distance between each pair of adjacent skeleton points in the order of the points in the minimum spanning tree of the main stem skeleton. Sum the distances between all pairs of adjacent skeleton points to obtain the sum of the distances as the main stem length H.
[0122] B. Ground diameter of Masson pine:
[0123] S11. Take the skeleton point at a preset fixed distance from the lowest skeleton point in the minimum spanning tree of the candidate points as the starting point, draw a line connecting the starting point and the skeleton point adjacent to the minimum spanning tree of the candidate points above it as the lower normal, and generate a plane perpendicular to the lower normal through the starting point as the lower plane.
[0124] The preset fixed distance is 0.2-0.3cm. S12. Take the nearest skeleton point along the minimum spanning tree of the candidate points above the starting point as the secondary starting point, and draw a line connecting the secondary starting point and the nearest skeleton point along the minimum spanning tree of the candidate points above it as the upper normal. Generate a plane perpendicular to the upper normal through the secondary starting point as the upper plane.
[0125] S13. In the main stem point cloud, all points in the main stem point cloud located between the upper and lower planes are vertically projected onto the lower plane to obtain projection points. The maximum distance between any two projection points is calculated as the ground diameter of the Masson pine.
[0126] C. Average blade length:
[0127] S21. Select a skeleton point from the minimum spanning tree of candidate points as the starting point, draw a line connecting the starting point and the skeleton point adjacent to the minimum spanning tree of candidate points above it as the lower normal, and generate a plane perpendicular to the lower normal through the starting point as the lower plane.
[0128] The preset fixed distance is 0.2-0.3cm.
[0129] S22. Take the nearest skeleton point of the minimum spanning tree of the candidate points above the starting point as the secondary starting point, and draw a line connecting the secondary starting point and the nearest skeleton point of the minimum spanning tree of the candidate points above it as the upper normal. Generate a plane perpendicular to the upper normal through the secondary starting point as the upper plane.
[0130] S23. In the main stem point cloud, all points in the main stem point cloud located between the upper and lower planes are projected vertically onto the lower plane to obtain projection points. A minimum α-bounding box is established to enclose all projection points. The average distance between all points on the boundary of the minimum α-bounding box and the projection point (only one point) projected onto the lower plane from the starting point is calculated:
[0131] If the average distance is greater than the preset distance threshold, the average distance will be used as the average leaf length of the region where the starting point is located on the Masson pine.
[0132] If the average distance is less than or equal to the preset distance threshold, it is assumed that there are no blades in the area where the starting point is located.
[0133] S24. Repeat steps S11-S14 for each skeleton point in the minimum spanning tree of the candidate points. Initially, take the lowest skeleton point in the minimum spanning tree of the candidate points as the starting point. After each traversal, shift the starting point upwards by one skeleton point along the minimum spanning tree to obtain the average leaf length of the region where each skeleton point is located.
[0134] D. Number of leaves:
[0135] S31. Select a skeleton point from the minimum spanning tree of candidate points as the starting point, draw a line connecting the starting point and the skeleton point adjacent to the minimum spanning tree of candidate points above it as the lower normal, and generate a plane perpendicular to the lower normal through the starting point as the lower plane.
[0136] The preset fixed distance is 0.2-0.3cm.
[0137] S32. Take the nearest skeleton point of the minimum spanning tree of the candidate points above the starting point as the secondary starting point, and draw a line connecting the secondary starting point and the nearest skeleton point of the minimum spanning tree of the candidate points above it as the upper normal. Generate a plane perpendicular to the upper normal through the secondary starting point as the upper plane.
[0138] S33. In the main stem point cloud, the number of points in the leaf point cloud located between the upper and lower planes.
[0139] If the number of points is greater than the preset number of points threshold, the average distance density of all points in the leaf point cloud located between the upper and lower planes is calculated, and the number of leaves in the starting area is obtained by multiplying the average distance density by a parameter.
[0140] If the number of points is not greater than the preset point threshold, then there are no leaves in the region where the starting point is located;
[0141] S34. Repeat steps S11-S14 for each skeleton point in the minimum spanning tree of the candidate points. Initially, take the lowest skeleton point in the minimum spanning tree of the candidate points as the starting point. After each traversal, shift the starting point upwards by one skeleton point along the minimum spanning tree to obtain the number of leaves in the region where each skeleton point is located.
[0142] The number of blades in the region of each skeleton point is calculated using the above method, and the results are summed to obtain the total number of blades.
[0143] This invention proposes a method for estimating the number of leaf blades: the density method. Existing point cloud density estimation methods include distance-based methods and block-based methods. Block-based methods divide the point cloud into blocks based on its distribution and count the number of points within each block to estimate the point cloud density. Distance-based average density representation estimates the sparseness of the point cloud distribution by calculating the average distance between points in the point cloud. However, since the specific distribution of the point cloud is unknown, it is difficult to select a block standard. Therefore, this experiment uses the distance-based density estimation method.
[0144] The specific method is as follows: The distance between points is represented as the distance between a point p in the point cloud and the nearest point in the point cloud to p. In a point cloud C with N points, let dis(p,q) represent the distance between p and any other point q, and dp represent the minimum distance between p and other points. Then:
[0145]
[0146] The average distance density of the leaf point cloud is:
[0147]
[0148] Where d represents the average distance density of the leaf point cloud;
[0149] As can be seen from the above formula, the smaller the average distance d, the denser the point cloud distribution; the larger d, the sparser the point cloud distribution, and the lower the density.
[0150] This experiment calculated the density of the point cloud after removing the main stem, leaving only the leaves. The point cloud density of the main stem affects the estimation of the leaf density. The average point cloud density of ten Masson pines with the main stem removed was calculated, and the number of leaves was manually counted. A certain linear relationship was found between the two. Therefore, by calculating the point cloud density of Masson pines with the main stem removed and multiplying it by a coefficient, the number of leaves can be estimated.
[0151] Experiments on four phenotypic parameters revealed that the accuracy rates for plant height, stem diameter (ground stem), and main stem length all exceeded 90%, meeting the project's eligibility criteria. Given the large number of leaves in Masson pine, the accuracy rate of estimation using the density method was acceptable.
[0152] In the embodiments of the present invention, those skilled in the art will also understand that all or part of the steps in the methods of the above embodiments can be implemented by a program instructing related hardware. The program can be stored in a computer-readable storage medium, including ROM / RAM, disk, optical disk, etc.
[0153] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.
Claims
1. An automatic measurement method for partial phenotypic parameters of Pinus massoniana based on three-dimensional point cloud images, characterized in that, The method includes the following steps: 1) Use a Kinect camera and a manual turntable to acquire three-dimensional point cloud images of Masson pine plants. During the acquisition process, rotate the Masson pine plants and record one three-dimensional point cloud image for each rotation. Perform N rotations on each Masson pine plant to record N three-dimensional point cloud images. 2) Preprocess each 3D point cloud image to obtain the 3D point cloud of the Masson pine plant itself; 3) Register the N preprocessed 3D point clouds to obtain a fused point cloud; 4) Extract the minimum spanning tree of candidate points and the skeleton points of the main stem, and then obtain the point cloud of the main stem and the point cloud of the leaves; Step 4) specifically involves: 4.1) The fused point cloud is sliced along the growth direction for the first time, and Euclidean clustering is performed on the point cloud in each slice layer obtained in the first slice to obtain each cluster. 4.2) Based on the clusters obtained in step 4.1), make the following judgments: For clusters in step 4.1) where the number of points is less than or equal to the preset threshold, the centroid of the cluster is directly extracted as candidate points for the main stem skeleton. For clusters with more than a preset threshold number of points in step 4.1), a second slice is made in any direction perpendicular to the growth direction. The points in each slice layer obtained in the second slice are clustered according to their gray values to obtain each cluster. The average gray value of each cluster is then calculated. The centroid of the slice layer of the second slice containing the cluster with the highest gray value is taken as the candidate point of the main stem skeleton. The set of candidate points is composed of all the candidate points of the main stem skeleton; 4.3) Construct the minimum spanning tree of the candidate point set, and then use the DAG single-source longest path algorithm to search for the longest path in the minimum spanning tree of the candidate points. All candidate points of the longest path form a set of main stem skeleton points containing the canopy. 4.4) Search downwards from the top candidate point in the set of main stem skeleton points containing the canopy. When a candidate point is found, if the search path length is greater than the canopy length threshold, stop the search. All candidate points found in the longest path up to the top candidate point are discarded. The remaining candidate points in the longest path are used as the skeleton points of the main stem. 4.5) Interpolate between the skeleton points of each main stem, inserting a point at fixed intervals; after completing the interpolation of the skeleton points of the main stem, use the points obtained by searching along the skeleton points from bottom to top using the radius search method of Kdtree as the points in the main stem point cloud, thereby obtaining the main stem point cloud; 4.6) Subtract the main stem point cloud from the fused point cloud obtained after registration in step 3) to obtain the leaf point cloud, thus achieving stem-leaf separation of Masson pine. 5) Obtain the main stem length H, diameter at ground level of Masson pine, average leaf length and number of leaves by processing the minimum spanning tree of candidate points, main stem point cloud and leaf point cloud.
2. The method for automatically measuring partial phenotypic parameters of Pinus massoniana based on three-dimensional point cloud images according to claim 1, characterized in that: The preprocessing in step 2) specifically involves: 2.1) Use a pass-through filter to remove the background from the 3D point cloud image; 2.2) The RANSAC plane fitting method is used to detect the desktop plane in the 3D point cloud image after removing the background, and then the desktop plane is removed using the plane segmentation method; 2.3) A radius-based filtering method is used to remove noise points from the 3D point cloud image after removing the background and the desktop plane; 2.4) Obtain the lowest point of the pot where the Masson pine plant is located from the 3D point cloud image after removing noise points. Then, estimate the soil plane height of the pot based on the lowest point and the pot height to determine the soil plane. Take the area above the soil plane as the region of interest and remove the area below the soil plane using a pass-through filter. All points in the region of interest constitute the 3D point cloud of the Masson pine plant itself.
3. The method for automatically measuring partial phenotypic parameters of Pinus massoniana based on three-dimensional point cloud images according to claim 1, characterized in that: Step 3) specifically refers to: 3.1) For both the source and target point clouds, traverse each point and perform the following processing to obtain the feature points of the source and target point clouds: Points in the point cloud are selected as processing points. A region point cloud is formed by considering all points within a neighborhood radius of the processing point. Principal component analysis (PCA) is used to estimate the covariance matrix of the region point cloud. The eigenvector corresponding to the smallest eigenvalue of the covariance matrix is then used as the normal vector of the processing point. Next, the angles between the normal vector of the processing point and the normal vectors of all points within the region point cloud and its neighbors are calculated. The mean of all these angles is then calculated as the average angle. An angle threshold is set, and the following judgments are made: If the mean angle of the processing point is greater than the angle threshold, the processing point is retained as a feature point; otherwise, the processing point is not retained as a feature point. 3.2) Calculate the feature descriptors of feature points in the source point cloud and the target point cloud, and use the feature descriptors in the source point cloud and the target point cloud to perform coarse registration using the sampling consistency method; 3.3) The source point cloud and target point after coarse registration coordinate transformation are used as input, and the nearest point iterative algorithm ICP method is used for fine registration. In the process of iteratively searching for the nearest point in the ICP method, the KD-tree based neighbor search algorithm is used for improvement to obtain the transformation matrix between the source point cloud and the target point cloud. 3.4) Take the last 3D point cloud in the N preprocessed 3D point clouds as the target point cloud, and the remaining 3D point clouds as the source point clouds. Traverse each source point cloud and repeat the above steps 3.1) to 3.4) to perform registration processing between the source point clouds and the target point clouds to obtain the transformation matrix between each pair of source point clouds and target point clouds. 3.5) Register all source point clouds to the target point cloud using their respective transformation matrices to obtain each registered source point cloud. Combine the registered source point clouds and target point clouds into a complete fused point cloud by point cloud addition.
4. The method for automatically measuring partial phenotypic parameters of Pinus massoniana based on three-dimensional point cloud images according to claim 1, characterized in that: In step 5), the length H of the main stem is obtained as follows: Based on the skeleton points of the main stem, a minimum spanning tree is built from the skeleton points of the main stem as the minimum spanning tree of the main stem skeleton. The distance between each two adjacent skeleton points is extracted according to the point order in the minimum spanning tree of the main stem skeleton. The sum of the distances between all two adjacent skeleton points is obtained as the length H of the main stem.
5. The method for automatically measuring partial phenotypic parameters of Pinus massoniana based on three-dimensional point cloud images according to claim 1, characterized in that: In step 5), the diameter at ground level of the Masson pine is obtained in the following manner: S11. Take the skeleton point at a preset fixed distance from the lowest skeleton point in the minimum spanning tree of the candidate points as the starting point, draw a line connecting the starting point and the skeleton point adjacent to the minimum spanning tree of the candidate points above it as the lower normal, and generate a plane perpendicular to the lower normal through the starting point as the lower plane. S12. Take the nearest skeleton point along the minimum spanning tree of the candidate points above the starting point as the secondary starting point, and draw a line connecting the secondary starting point and the nearest skeleton point along the minimum spanning tree of the candidate points above it as the upper normal. Generate a plane perpendicular to the upper normal through the secondary starting point as the upper plane. S13. Project all points in the main stem point cloud located between the upper and lower planes vertically to the lower plane to obtain projection points, and calculate the maximum distance between every two projection points as the ground diameter of the Masson pine.
6. The method for automatically measuring partial phenotypic parameters of Pinus massoniana based on three-dimensional point cloud images according to claim 1, characterized in that: In step 5), the average blade length is obtained in the following manner: S21. Select a skeleton point from the minimum spanning tree of candidate points as the starting point, draw a line connecting the starting point and the skeleton point adjacent to the minimum spanning tree of candidate points above it as the lower normal, and generate a plane perpendicular to the lower normal through the starting point as the lower plane. S22. Take the nearest skeleton point of the minimum spanning tree of the candidate points above the starting point as the secondary starting point, and draw a line connecting the secondary starting point and the nearest skeleton point of the minimum spanning tree of the candidate points above it as the upper normal. Generate a plane perpendicular to the upper normal through the secondary starting point as the upper plane. S23. Project all points in the main stem point cloud located between the upper and lower planes vertically to the lower plane to obtain projection points. Construct the minimum α-bounding box that encloses all projection points. Calculate the average distance between all points on the boundary of the minimum α-bounding box and the projection points projected onto the lower plane from the starting point: If the average distance is greater than the preset distance threshold, the average distance is used as the average leaf length of the region where the starting point is located. If the average distance is less than or equal to the preset distance threshold, it is assumed that there are no blades in the area where the starting point is located. S24. Repeat steps S11-S14 for each skeleton point in the minimum spanning tree of the candidate points to obtain the average leaf length of the region where each skeleton point is located.
7. The method for automatically measuring partial phenotypic parameters of Pinus massoniana based on three-dimensional point cloud images according to claim 1, characterized in that: In step 5), the number of blades is obtained in the following manner: S31. Select a skeleton point from the minimum spanning tree of candidate points as the starting point, draw a line connecting the starting point and the skeleton point adjacent to the minimum spanning tree of candidate points above it as the lower normal, and generate a plane perpendicular to the lower normal through the starting point as the lower plane. S32. Take the nearest skeleton point of the minimum spanning tree of the candidate points above the starting point as the secondary starting point, and draw a line connecting the secondary starting point and the nearest skeleton point of the minimum spanning tree of the candidate points above it as the upper normal. Generate a plane perpendicular to the upper normal through the secondary starting point as the upper plane. S33. Count the number of points in the leaf point cloud located between the upper and lower planes; If the number of points is greater than the preset number of points threshold, the average distance density of all points in the leaf point cloud located between the upper and lower planes is calculated, and the number of leaves in the starting area is obtained by multiplying the average distance density by a parameter. If the number of points is not greater than the preset point threshold, then there are no leaves in the region where the starting point is located; S34. Repeat steps S11-S14 for each skeleton point in the minimum spanning tree of the candidate points to obtain the number of leaves in the region of each skeleton point.