CN109697731B - Plant canopy volume calculation method and system based on point cloud data and related equipment - Google Patents
Plant canopy volume calculation method and system based on point cloud data and related equipment Download PDFInfo
- Publication number
- CN109697731B CN109697731B CN201811571128.8A CN201811571128A CN109697731B CN 109697731 B CN109697731 B CN 109697731B CN 201811571128 A CN201811571128 A CN 201811571128A CN 109697731 B CN109697731 B CN 109697731B
- Authority
- CN
- China
- Prior art keywords
- cross
- section
- plant canopy
- calculating
- points
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 64
- 230000011218 segmentation Effects 0.000 claims abstract description 34
- 230000006870 function Effects 0.000 claims description 38
- 238000000034 method Methods 0.000 claims description 31
- 238000012545 processing Methods 0.000 claims description 4
- 238000010586 diagram Methods 0.000 description 18
- 230000000694 effects Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 241000283973 Oryctolagus cuniculus Species 0.000 description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 229910052799 carbon Inorganic materials 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10028—Range image; Depth image; 3D point clouds
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30188—Vegetation; Agriculture
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Image Generation (AREA)
Abstract
The invention relates to the technical field of volume calculation, in particular to a plant canopy volume calculation method and system based on point cloud data and related equipment, and aims to solve the problem that the plant canopy volume calculation in the prior art is not ideal. The calculation method of the invention comprises the following steps: projecting the coordinate information of the plant canopy point cloud data to a z axis by taking the three-dimensional coordinate data of the plant canopy point cloud data as an initial value, and calculating a z coordinate probability density function; calculating a curvature local extreme point of the z coordinate probability density function according to the z coordinate probability density function; for each local extreme point, drawing a plane as a cross section through the point and perpendicular to the z axis; according to the preset number of segmentation layers, segmenting among the local extreme points to obtain a plurality of cross sections perpendicular to the z axis; and dividing the plant canopy into a cone and a plurality of platforms from top to bottom according to the cross section, and further calculating the volume of the plant canopy. The invention improves the accuracy of the volume calculation of the plant canopy.
Description
Technical Field
The invention relates to the technical field of volume calculation, in particular to a plant canopy volume calculation method and system based on point cloud data and related equipment.
Background
Accurate measurement of plant canopy volume has long been a problem for researchers due to the irregularity and complexity of the plant canopy spatial structure itself. The traditional method is to divide the plant canopy at equal intervals and abstract the plant canopy into a single geometric body, and the method has great limitation and great error.
Most of experts and scholars at home and abroad adopt an equal-interval slicing method to abstract a cross section into two-dimensional geometric figures such as a circle, an ellipse, a fan, a trapezoid and the like to calculate the volume of a crown, but the calculation of the cross section and the calculation of a side area by the research results cannot ensure high calculation precision.
Disclosure of Invention
In order to solve the problems in the prior art, the invention provides a plant canopy volume calculation method and system based on point cloud data and related equipment, and the precision of the plant canopy volume calculation is improved.
In a first aspect of the present invention, a plant canopy volume calculation method based on point cloud data is provided, the method includes:
projecting the coordinate information of the plant canopy point cloud data to a z axis by taking the three-dimensional coordinate data of the plant canopy point cloud data as an initial value, and calculating a z coordinate probability density function;
calculating a curvature local extreme point of the z coordinate probability density function according to the z coordinate probability density function;
for each local extreme point, drawing a plane as a cross section through the point and perpendicular to the z-axis; according to a preset number of segmentation layers, segmenting among the local extreme points to obtain a plurality of cross sections perpendicular to the z axis;
and dividing the plant canopy into a cone and a plurality of platforms from top to bottom according to the cross section, and further calculating the volume of the plant canopy.
Preferably, the step of calculating the curvature local extreme point of the z-coordinate probability density function according to the z-coordinate probability density function includes:
the probability density is corrected as follows:
p'＝corrp*p
wherein,
max (z), min (z) respectively represent the maximum value and the minimum value of the projection of the plant canopy point cloud data on the z axis, and max (p) represents the maximum value of the probability density p;
the curvature of the discrete points on the z-axis projection is calculated as follows:
wherein,
△fk(t)＝f(t+k△t)-f(t+(k-1)△t)
c (t) is the calculated curvature value, f (t) is a function, k is 1,2,3, Δ t is the unit step length of the parameter t of the function f (t), the symbol "x" represents the cross product, the symbol "·" represents the dot product, and the symbol "| |" represents the modular length;
selecting a point with a curvature value larger than a curvature threshold value from discrete points on the z-axis projection as an undetermined extreme value point;
and selecting the point with the maximum curvature value from the adjacent undetermined extreme points to obtain a local extreme point.
Preferably, the step of dividing between the local extreme points according to a preset number of dividing layers to obtain a plurality of cross sections perpendicular to the z-axis includes:
calculating the segmentation step length between each pair of adjacent local extreme points according to the following formula:
performing equal-distance segmentation between each pair of adjacent local extreme points according to the corresponding segmentation step length to obtain a plurality of cross sections perpendicular to the z axis;
wherein, LEziAnd LEzi+1The heights of the ith and (i + 1) th local extreme points are respectively the heights of the ith and (i + 1) th local extreme points, the ith and (i + 1) th local extreme points are two adjacent local extreme points, i is 1,2, …, m, and m is the number of the local extreme points; and delta z is the height of the point cloud, and L is the preset segmentation layer number.
Preferably, the step of dividing the plant canopy into a vertebral body and a plurality of table bodies from top to bottom according to the cross section and further calculating the volume of the plant canopy comprises:
respectively projecting points in the upper and lower specific distance ranges of each cross section onto the corresponding cross section, and further calculating boundary points of each cross section;
connecting the boundary points of each cross section in sequence to obtain a polygonal contour of the corresponding cross section, and further calculating the area of each cross section;
dividing the plant canopy into a cone and a plurality of platforms from top to bottom by taking the cross section as an interface;
and respectively calculating the volumes of the vertebral body and the platform body according to the area of the cross section, thereby obtaining the volume of the plant canopy.
Preferably, the step of "projecting points in a specific distance range above and below each cross section onto the corresponding cross section respectively, and calculating boundary points of each cross section" includes:
calculating the average distance between points of the plant canopy point cloud data according to the following formula:
wherein,
rho is the density of the point cloud in the height direction, N is the number of the points in the point cloud data, and S is the projection area of the point cloud in the horizontal direction of the geodetic coordinate; delta z is the point cloud height;
respectively projecting points in a specific distance range above and below each cross section onto the corresponding cross section, wherein the specific distance isz/2；
And calculating boundary points of the cross section by adopting an alpha-shape algorithm of a two-dimensional plane boundary contour for the point cloud data on each cross section.
Preferably, the step of "sequentially connecting boundary points of each cross section respectively to obtain a polygonal contour of the corresponding cross section, and then calculating an area of each cross section" includes:
connecting the boundary points of each cross section in sequence to obtain a polygonal profile of the corresponding cross section;
the first vertex p of the polygon outline1Respectively connecting with other vertexes to obtain a plurality of triangles;
calculating the area of the polygonal contour as the area of the corresponding cross section according to the following formula:
wherein S isiIs a triangle p1pipi+1R is the number of vertices of the polygonal profile; a. theiIs a sign coefficient, if p1、pi、pi+1In a counterclockwise arrangement, then AiIf 1, if p1、pi、pi+1In clockwise arrangement, Ai＝-1。
Preferably, the step of calculating the volume of the vertebral body and the volume of the table body respectively according to the area of the cross section so as to obtain the volume of the plant canopy comprises the following steps:
the volume of the ith table was calculated as follows:
the volume of the vertebral body is calculated as follows:
the volume of the plant canopy was calculated as follows:
wherein S isiAnd Si+1The areas of the cross sections corresponding to the lower bottom surface and the upper bottom surface of the ith platform body are respectively, i is the serial number of the platform body and is arranged from the bottom of the plant canopy to the top in sequence, i is 1,2iAnd zi+1Respectively is the z-axis coordinate of the lower bottom surface and the upper bottom surface of the ith platform body; sg-1Is the area of the cross section corresponding to the bottom surface of the vertebral body, zgAnd zg-1Respectively the z-axis coordinates of the base and apex of the vertebral body.
In a second aspect of the present invention, a plant canopy volume calculation system based on point cloud data is provided, the system comprising:
the probability density function calculation module is used for projecting the coordinate information of the plant canopy point cloud data to a z axis by taking the three-dimensional coordinate data of the plant canopy point cloud data as an initial value, and calculating a z coordinate probability density function;
the extreme point calculating module is used for calculating a curvature local extreme point of the z coordinate probability density function according to the z coordinate probability density function;
the segmentation module is used for drawing a plane as a cross section for each local extreme point through the point and perpendicular to the z axis; according to a preset number of segmentation layers, segmenting among the local extreme points to obtain a plurality of cross sections perpendicular to the z axis;
and the volume calculation module is used for dividing the plant canopy into a cone and a plurality of platforms from top to bottom according to the cross section, and further calculating the volume of the plant canopy.
In a third aspect of the invention, a storage device is proposed, comprising a memory adapted to store a program adapted to be loaded and executed by a processor to implement the plant canopy volume calculation method based on point cloud data as described above.
In a fourth aspect of the present invention, a processing apparatus is provided, comprising:
a processor adapted to load a program;
a memory adapted to store the program;
the program is adapted to be loaded and executed by the processor to implement the plant canopy volume calculation method based on point cloud data as described above.
Compared with the closest prior art, the invention has the following beneficial effects:
according to the plant canopy volume calculation method based on the point cloud data, the density information of the three-dimensional point cloud in the vertical direction is utilized, an unequal-interval segmentation strategy is provided, meanwhile, the alpha-shape algorithm of solving the boundary contour by using the two-dimensional plane is utilized, a more accurate cross section contour is provided, and the defect that the cross section contour and the side shape calculation of the traditional equidistant slicing method are rough is overcome. The method is applied to the calculation of the volume of the plant canopy, then the plant canopy is divided into a combined body of a cone at the top and a plurality of platforms at the bottom, and finally the volume of the plant canopy is calculated by adopting an accumulative method of the volumes of the platforms and the volumes of the cones, so that the accuracy of the calculation of the volume of the plant canopy is improved. The plant canopy volume calculation method can be widely applied to the research of the agricultural and forestry fields, and can also be applied to the environmental evaluation of urban greening conditions, carbon emission and the like.
Drawings
FIGS. 1(a) and 1(b) are schematic diagrams of point cloud data of a rabbit and a crown, respectively;
FIG. 2 is a schematic diagram of the main steps of the embodiment of the plant canopy volume calculation method based on point cloud data;
FIG. 3 is a schematic diagram illustrating a probability density function of point cloud data of a crown according to an embodiment of the present invention;
FIG. 4 is a schematic diagram comparing the effect of the prior art segmentation method and the segmentation method of the present invention;
FIG. 5 is a schematic diagram illustrating comparison of cross-sectional profiles of an embodiment of the present invention using an alpha-shape algorithm with four other strategies;
FIG. 6 is a diagram illustrating an embodiment of the present invention in which a polygon is divided into a plurality of triangles;
FIG. 7 is a schematic diagram of the plant canopy divided into a cone and a plurality of platforms from top to bottom according to the embodiment of the present invention;
FIG. 8 is a schematic diagram of the main components of the plant canopy volume calculation system based on point cloud data according to the embodiment of the present invention.
Detailed Description
Preferred embodiments of the present invention are described below with reference to the accompanying drawings. It should be understood by those skilled in the art that these embodiments are only for explaining the technical principle of the present invention, and are not intended to limit the scope of the present invention.
Fig. 1(a) and 1(b) are schematic diagrams of point cloud data of a rabbit and a crown, respectively. The point cloud data shown in fig. 1(a) is data of the rabbit body surface obtained by laser scanning. The point cloud data shown in fig. 1(b) is obtained by scanning a laser beam and is a crown point cloud data. Because of the gaps among the leaves, when the tree leaves are irradiated by light, a part of the light can penetrate through the plant canopy, and a part of the light can be blocked by the branches or the leaves and returns after being irradiated into the plant canopy. Therefore, when the plant canopy is scanned by the instrument, the obtained point cloud includes the data inside the canopy in addition to the surface data. Therefore, for the particularity of the plant canopy, the invention adopts an unequal interval segmentation method based on density and uses an alpha-shape algorithm to detect the contour line of the cross section so as to more accurately calculate the volume of the plant canopy.
FIG. 2 is a schematic diagram of the main steps of the embodiment of the plant canopy volume calculation method based on point cloud data. As shown in FIG. 2, the calculation method of the present embodiment includes steps S1-S4:
and step S1, projecting the coordinate information of the plant canopy point cloud data to a z-axis by taking the three-dimensional coordinate data of the plant canopy point cloud data as an initial value, and calculating a z-coordinate probability density function.
FIG. 3 is a schematic diagram of a probability density function of point cloud data of a crown according to an embodiment of the present invention. As shown in fig. 3, the left half is a schematic diagram of a crown, and the right half is a probability density function corresponding to the point cloud data of the crown. The vertical axis is the z coordinate of the point cloud data, and the horizontal axis is the probability density of the point cloud data projected to the z axis.
And step S2, calculating the curvature local extreme point of the z coordinate probability density function according to the z coordinate probability density function. The method may specifically include steps S21-S24:
step S21, the probability density is corrected according to the formulas (1) to (2):
p'＝corrp*p(1)
wherein, max (z), min (z) respectively represent the maximum value and the minimum value of the plant canopy point cloud data projected on the z-axis, and max (p) represents the maximum value of the probability density p.
Step S22, calculating the curvature of the discrete point on the z-axis projection according to the formulas (3) to (4):
△fk(t)＝f(t+k△t)-f(t+(k-1)△t) (4)
where c (t) is the calculated curvature value, f (t) is a function, k is 1,2,3, Δ t is the unit step length of the parameter t of the function f (t), the symbol "x" represents a cross product, the symbol "·" represents a dot product, and the symbol "| |" represents a modular length.
The idea of calculating the curvature of discrete points using differences is proposed in the literature [ Blankenburg C, Daul C, Ohser J.parameter free standardization of currents in 3D images [ C ]// IEEE International Conference on image processing. IEEE,2016: 1081-.
And step S23, selecting points with curvature values larger than a curvature threshold value from the discrete points on the z-axis projection as undetermined extreme value points.
When the selected curvature threshold is smaller, the number of local curvature extreme points is increased, so that the time cost is increased; when the curvature threshold is large, the number of local curvature extreme points is reduced, and the calculation accuracy is insufficient. According to the invention, through a large number of experiments, the corrected z coordinate probability density function graph can obtain a better local extreme point selection effect when the curvature threshold value is set to be 0.03.
And step S24, selecting the point with the maximum curvature value from the adjacent extreme points to be determined to obtain a local extreme point.
Step S3, drawing a plane as a cross section for each local extreme point through the point and perpendicular to the z axis; and according to the preset number of segmentation layers, segmenting among the local extreme points to obtain a plurality of cross sections perpendicular to the z axis. The method may specifically include steps S31-S33:
in step S31, for each local extreme point, a plane is drawn through the point and perpendicular to the z-axis as a cross-section.
Step S32, calculating a dividing step between each pair of adjacent local extreme points according to equation (5):
wherein, LEziAnd LEzi+1The heights of the ith and (i + 1) th local extreme points are respectively, the ith and (i + 1) th local extreme points are two adjacent local extreme points and are serial numbers of the local extreme points, and i is 1,2, …, m is sequentially arranged from bottom to top and is the number of the local extreme points; Δ z is the height of the point cloud, and L is the preset number of segmentation layers (20 in this example);
indicating rounding up.
And step S33, performing equal-distance segmentation between each pair of adjacent local extreme points according to the corresponding segmentation step size, thereby obtaining a plurality of cross sections perpendicular to the z axis.
According to the above algorithm, the number of segmentation layers obtained is usually larger than L for rounding up.
Fig. 4 is a schematic diagram comparing a prior art segmentation method with the segmentation method of the present invention. As shown in fig. 4, the left graph is a schematic diagram of equal-pitch segmentation in the prior art, and the right graph is a schematic diagram of unequal-pitch segmentation in the method. Wherein z is1、z2、z3、z4、z5Each representing a different cross-section. On the longitudinal section view, the intersection points of the cross section and the original contour are connected in sequence to obtain the divided contour, and the contour divided by the method is closer to the original contour and has smaller error.
And step S4, dividing the plant canopy into a cone and a plurality of platforms from top to bottom according to the cross section, and further calculating the volume of the plant canopy. The method may specifically include steps S41-S44:
and step S41, respectively projecting points in the upper and lower specific distance ranges of each cross section onto the corresponding cross section, and further calculating boundary points of each cross section. Steps S411 to S413 may be specifically included:
step S411, calculating the average distance between points of the plant canopy point cloud data according to the formulas (6) to (7):
wherein rho is the density of the point cloud in the height direction, N is the number of the points in the point cloud data, and S is the projection area of the point cloud in the horizontal direction of the geodetic coordinate; and delta z is the height of the point cloud.
Step S412, respectively projecting points in the upper and lower specific distance ranges of each cross section onto the corresponding cross section, wherein the specific distance isz/2。
Step S413, calculating boundary points of the cross section by using an alpha-shape algorithm of a two-dimensional plane boundary outline for the point cloud data on each cross section, and orderly storing the boundary points into a boundary point set.
And step S42, sequentially connecting the boundary points of each cross section respectively to obtain the polygonal contour of the corresponding cross section, and further calculating the area of each cross section. Steps S421 to S423 may be specifically included:
step S421, sequentially connecting the boundary points of each cross section to obtain a polygonal profile of the corresponding cross section.
FIG. 5 is a schematic diagram showing the comparison of the effect of the cross-sectional profile of the present embodiment using the alpha-shape algorithm with the cross-sectional profiles of the other four strategies. As shown in fig. 5, the first is a document "wangqi, huohong, wu yanlan, etc.. in a crown volume automatic calculation method [ J ] journal data based on point cloud data, academy of northwest forest, 2017,32(2):242 + 246 ], a slicing method is adopted to divide the crown point cloud data according to the z-axis, abstract the cross section into a circle, abstract each layer into a circular table or a cone, and finally, the volume is accumulated and solved; the second is a document "yellow key, Zhao \31066; xi, Shulong, etc.. laser measurement method and experimental study of fruit tree canopy volume [ C ]. Chinese agriculture engineering society 2011 academic annual meeting 2011", in which the cross section is abstracted to an ellipse; the third is a document 'three-dimensional green quantity estimation and analysis research in six regions of Tangxuehai, Beijing City [ D ]. Beijing university of forestry, 2011'; and the fourth method is a document 'Wangqi, Huhong, Wuyanlan, and the like', an automatic calculation method for the volume of a crown based on point cloud data [ J ]. proceedings of the northwest forest academy of academic, 2017,32(2):242 and 246 ], wherein the cross section is divided into a plurality of small sectors, and the area of the cross section is calculated by adopting a sector area approximation method. By contrast, the alpha-shape algorithm applied herein has the highest precision.
Step S422, the first vertex p of the polygon outline is processed1And connecting with other vertexes respectively to obtain a plurality of triangles.
Fig. 6 is a schematic diagram of the present embodiment in which a polygon is divided into a plurality of triangles. As shown in FIG. 6, the polygon has a total of 6 vertices, each being p1、p2、...、p6The first vertex p1Respectively connected with other vertexes to obtain four triangles p1p2p3、p1p3p4、p1p4p5、p1p5p6The areas of the four triangles are respectively marked as S2、S3、S4、S5。
Step S423, calculating the area of the polygonal contour as the area of the corresponding cross section according to the formula (8):
wherein S isiIs a triangle p1pipi+1R is the number of vertexes of the polygonal contour, namely the number of boundary points of the cross section corresponding to the polygonal contour; a. theiIs a sign coefficient, if p1、pi、pi+1Arranged counterclockwise on the polygonal outline (triangle p in fig. 6)1p2p3、p1p3p4、p1p4p5) Then A isiIf 1, if p1、pi、pi+1Arranged clockwise on the polygonal outline (triangle p in fig. 6)1p5p6) Then A isi＝-1。
And step S43, dividing the plant canopy into a cone and a plurality of platforms from top to bottom by taking the cross section as an interface.
And step S44, respectively calculating the volumes of the cone and the table according to the area of the cross section, and further obtaining the volume of the plant canopy. The method may specifically include steps S441-S443:
step S441, the volume of the i-th stage body is calculated according to the formula (9):
step S442, calculating the volume of the vertebral body according to the formula (10):
step S443, calculating the volume of the plant canopy according to equation (11):
wherein，SiAnd Si+1The cross section areas corresponding to the lower bottom surface and the upper bottom surface of the ith platform body (namely the lower bottom area and the upper bottom area of the ith platform body) are respectively, i is the serial number of the platform body and is arranged from the bottom of the plant canopy upwards, i is 1,2, g-2, and g is the number of the cross sections (including the cross section drawn at the uppermost extreme point, but the cross section has no meaning to the volume calculation, so the sum of the volume of g-2 platform bodies and one cone body is actually calculated), and z is the sum of the volume of the g-2 platform bodies and the volume of the cone body), and the cross section area is the sum of the volume of the plantiAnd zi+1Respectively is the z-axis coordinate of the lower bottom surface and the upper bottom surface of the ith platform body; sg-1Is the area of the cross section corresponding to the bottom surface of the vertebral body (i.e. the bottom area of the vertebral body), zgAnd zg-1Respectively, the z-axis coordinates of the base and apex of the vertebral body.
Fig. 7 is a schematic view of the plant canopy divided into a vertebral body and a plurality of platforms from top to bottom in this embodiment. As shown in FIG. 7, S1And S2The areas of the lower bottom surface and the upper bottom surface of the 1 st table body are respectively, namely the areas of the cross sections corresponding to the lower bottom surface and the upper bottom surface of the 1 st table body; z is a radical of1、z2Coordinates of the lower bottom surface and the upper bottom surface of the 1 st platform body on the z axis are respectively; in the same way, S2And S3Respectively the lower bottom area and the upper bottom area of the 2 nd stage body, Sg-1Is the area of the cross section corresponding to the bottom surface of the vertebral body, zgAnd zg-1Respectively, the z-axis coordinates of the base and apex of the vertebral body.
Compared with documents 'Wangqi, Huhong, Wuyanlan, and the like', the method for calculating the volume of the crown by unequal-interval segmentation based on the density extreme points automatically calculates the volume of the crown based on point cloud data [ J ]. the academy of northwest forest academy, 2017,32(2):242 + 246 ], and improves the calculation accuracy in the documents 'estimation and analysis research of the three-dimensional green volume in six areas of Tangxuehai, Beijing City, and [ D ]. Beijing forestry university, 2011'.
Although the foregoing embodiments describe the steps in the above sequential order, those skilled in the art will understand that, in order to achieve the effect of the present embodiments, the steps may not be executed in such an order, and may be executed simultaneously (in parallel) or in an inverse order, and these simple variations are within the scope of the present invention.
Based on the same technical concept as the embodiment of the calculation method, the invention also provides a calculation system, which is specifically described below.
FIG. 8 is a schematic diagram of the main components of the plant canopy volume calculation system based on point cloud data according to the embodiment of the present invention. As shown in fig. 8, the plant canopy volume calculation system 1 of the present embodiment includes: a probability density function calculation module 10, an extreme point calculation module 20, a segmentation module 30 and a volume calculation module 40.
The probability density function calculation module 10 is configured to project coordinate information of plant canopy point cloud data to a z-axis with three-dimensional coordinate data of the plant canopy point cloud data as an initial value, and calculate a z-coordinate probability density function; the extreme point calculating module 20 is configured to calculate a curvature local extreme point of the z-coordinate probability density function according to the z-coordinate probability density function; the segmentation module 30 is configured to draw a plane as a cross-section for each local extreme point passing through the point and perpendicular to the z-axis; according to the preset number of segmentation layers, segmenting among the local extreme points to obtain a plurality of cross sections perpendicular to the z axis; the volume calculation module 40 is used for dividing the plant canopy into a cone and a plurality of platforms from top to bottom according to the cross section, and further calculating the volume of the plant canopy.
In this embodiment, the extreme point calculating module 20 may specifically include: the device comprises a correction unit, a curvature calculation unit, an undetermined extreme point selection unit and a local extreme point selection unit.
Wherein, the correcting unit is used for correcting the probability density according to the formulas (1) to (2); the curvature calculation unit is used for calculating the curvature of the discrete point on the z-axis projection according to the formulas (3) to (4); the undetermined extreme point selecting unit is used for selecting a point with a curvature value larger than a curvature threshold value from discrete points on the z-axis projection as an undetermined extreme point; the local extreme point selecting unit is used for selecting a point with the maximum curvature value from the adjacent undetermined extreme points to obtain a local extreme point.
In this embodiment, the dividing module 30 may specifically include: the device comprises an extreme point cross section generating unit, a segmentation step calculating unit and a segmentation unit.
The extreme point cross section generating unit is used for drawing a plane as a cross section for each local extreme point, passing through the point and being perpendicular to the z axis; the dividing step calculating unit is used for calculating the dividing step between each pair of adjacent local extreme points according to the formula (5); the dividing unit is used for performing equal-interval division between each pair of adjacent local extreme points according to the corresponding dividing step length, and therefore a plurality of cross sections perpendicular to the z axis are obtained.
In this embodiment, the volume calculating module 40 may specifically include: the device comprises a boundary point calculation unit, a cross section area calculation unit, a dividing unit and a volume calculation unit.
The boundary point calculation unit is used for projecting points in the upper specific distance range and the lower specific distance range of each cross section onto the corresponding cross section respectively so as to calculate the boundary point of each cross section; the cross section area calculation unit is used for respectively and sequentially connecting the boundary points of each cross section to obtain a polygonal contour of the corresponding cross section, and further calculating the area of each cross section; the dividing unit is used for dividing the plant canopy into a cone and a plurality of platforms from top to bottom by taking the cross section as an interface; the volume calculation unit is used for respectively calculating the volumes of the vertebral body and the frustum body according to the area of the cross section, and further obtaining the volume of the plant canopy.
Specifically, the boundary point calculating unit in this embodiment may include: calculation subunit for average distance between points, projection subunit, and boundary point calculation subunit
The inter-point average distance calculation subunit is used for calculating the inter-point average distance of the plant canopy point cloud data according to the formulas (6) to (7); the projection subunit is used for projecting the points in the upper and lower specific distance ranges of each cross section onto the corresponding cross section respectively; and the boundary point calculation subunit is used for calculating the boundary points of the cross section by adopting an alpha-shape algorithm of a two-dimensional plane boundary contour for the point cloud data on each cross section.
Specifically, the cross-sectional area calculating unit in this embodiment may include: the device comprises a polygon outline generating subunit, a triangle generating subunit and a cross-sectional area calculating subunit.
The polygon outline generating subunit is used for respectively and sequentially connecting the boundary points of each cross section to obtain the polygon outline of the corresponding cross section; the triangle generation subunit is used for generating the first vertex p of the polygon outline1Respectively connecting with other vertexes to obtain a plurality of triangles; the cross-sectional area calculating subunit is configured to calculate an area of the polygonal contour as an area of the corresponding cross-section according to equation (8).
Specifically, the volume calculating unit in this embodiment may include: a table volume calculating subunit, a vertebral volume calculating subunit and a summing subunit.
The table body volume calculating subunit is used for calculating the volume of the table body according to a formula (9); the vertebral body volume calculating subunit is used for calculating the volume of the vertebral body according to the formula (10); the summation subunit is used to calculate the volume of the plant canopy according to equation (11).
In the present application, the computing system is divided into modules and units only for better understanding of the functions related to the technical solutions of the present invention, and in practice, the functions corresponding to the modules may be loaded and executed by a single or multiple hardware.
Based on the above calculation method, the present invention further proposes an embodiment of a storage device, comprising a memory adapted to store a program adapted to be loaded and executed by a processor to implement the above-mentioned plant canopy volume calculation method based on point cloud data.
Further, the present invention also provides a processing apparatus, comprising: a processor and a memory. Wherein the processor is adapted to load a program, the memory is adapted to store said program, said program is adapted to be loaded and executed by said processor to implement the plant canopy volume calculation method based on point cloud data as described above.
Those of skill in the art will appreciate that the method steps of the examples described in connection with the embodiments disclosed herein may be embodied in electronic hardware, computer software, or combinations of both, and that the components and steps of the examples have been described above generally in terms of their functionality in order to clearly illustrate the interchangeability of electronic hardware and software. Whether such functionality is implemented as electronic hardware or software depends upon the particular application and design constraints imposed on the solution. Skilled artisans may implement the described functionality in varying ways for each particular application, but such implementation decisions should not be interpreted as causing a departure from the scope of the present invention.
So far, the technical solutions of the present invention have been described in connection with the preferred embodiments shown in the drawings, but it is easily understood by those skilled in the art that the scope of the present invention is obviously not limited to these specific embodiments. Equivalent changes or substitutions of related technical features can be made by those skilled in the art without departing from the principle of the invention, and the technical scheme after the changes or substitutions can fall into the protection scope of the invention.
Claims (9)
1. A plant canopy volume calculation method based on point cloud data is characterized by comprising the following steps:
projecting the coordinate information of the plant canopy point cloud data to a z axis by taking the three-dimensional coordinate data of the plant canopy point cloud data as an initial value, and calculating a z coordinate probability density function;
calculating a curvature local extreme point of the z coordinate probability density function according to the z coordinate probability density function;
for each local extreme point, drawing a plane as a cross section through the point and perpendicular to the z-axis; according to a preset number of segmentation layers, segmenting among the local extreme points to obtain a plurality of cross sections perpendicular to the z axis;
dividing the plant canopy into a cone and a plurality of platforms from top to bottom according to the cross section, and further calculating the volume of the plant canopy;
the method for calculating the volume of the plant canopy comprises the following steps: respectively projecting points in the upper and lower specific distance ranges of each cross section onto the corresponding cross section, and further calculating boundary points of each cross section; connecting the boundary points of each cross section in sequence to obtain a polygonal contour of the corresponding cross section, and further calculating the area of each cross section; dividing the plant canopy into a cone and a plurality of platforms from top to bottom by taking the cross section as an interface; and respectively calculating the volumes of the vertebral body and the platform body according to the area of the cross section, thereby obtaining the volume of the plant canopy.
2. The method of claim 1, wherein the step of calculating the curvature local extreme point of the z-coordinate probability density function according to the z-coordinate probability density function comprises:
the probability density is corrected as follows:
p'＝corrp*p
wherein,
max (z), min (z) respectively represent the maximum value and the minimum value of the projection of the plant canopy point cloud data on the z axis, and max (p) represents the maximum value of the probability density p;
the curvature of the discrete points on the z-axis projection is calculated as follows:
wherein,
Δfk(t)＝f(t+kΔt)-f(t+(k-1)Δt)
c (t) is the calculated curvature value, f (t) is a function, k is 1,2,3, Δ t is the unit step length of the parameter t of the function f (t), the symbol "x" represents the cross product
The expression dot product, the symbol "| | |" represents the modular length;
selecting a point with a curvature value larger than a curvature threshold value from discrete points on the z-axis projection as an undetermined extreme value point;
and selecting the point with the maximum curvature value from the adjacent undetermined extreme points to obtain a local extreme point.
3. The method for calculating plant canopy volume according to claim 1, wherein the step of dividing between the local extreme points according to a preset number of dividing layers to obtain a plurality of cross sections perpendicular to the z-axis comprises:
calculating the segmentation step length between each pair of adjacent local extreme points according to the following formula:
performing equal-distance segmentation between each pair of adjacent local extreme points according to the corresponding segmentation step length to obtain a plurality of cross sections perpendicular to the z axis;
wherein, LEziAnd LEzi+1The heights of the ith and (i + 1) th local extreme points are respectively the heights of the ith and (i + 1) th local extreme points, the ith and (i + 1) th local extreme points are two adjacent local extreme points, i is 1,2, …, m, and m is the number of the local extreme points; and delta z is the height of the point cloud, and L is the preset number of segmentation layers.
4. The method of claim 1, wherein the step of projecting points in a specific distance range above and below each cross section onto the corresponding cross section to calculate the boundary points of each cross section comprises:
calculating the average distance between points of the plant canopy point cloud data according to the following formula:
wherein,
rho is the density of the point cloud in the height direction, N is the number of the points in the point cloud data, and S is the projection area of the point cloud in the horizontal direction of the geodetic coordinate; Δ z is the point cloud height;
respectively projecting points in a specific distance range above and below each cross section onto the corresponding cross section, wherein the specific distance isz/2；
And calculating boundary points of the cross section by adopting an alpha-shape algorithm of a two-dimensional plane boundary contour for the point cloud data on each cross section.
5. The method of claim 1, wherein the step of sequentially connecting the boundary points of each cross section to obtain the polygonal contour of the corresponding cross section and calculating the area of each cross section comprises:
connecting the boundary points of each cross section in sequence to obtain a polygonal profile of the corresponding cross section;
the first vertex p of the polygon outline1Respectively connecting with other vertexes to obtain a plurality of triangles;
calculating the area of the polygonal contour as the area of the corresponding cross section according to the following formula:
wherein S isiIs a triangle p1pipi+1R is the number of vertices of the polygonal profile; a. theiIs a sign coefficient, if p1、pi、pi+1In a counterclockwise arrangement, then AiIf 1, if p1、pi、pi+1In clockwise arrangement, Ai＝-1。
6. The method as claimed in claim 1, wherein the step of calculating the volume of the cone and the table according to the area of the cross section and obtaining the volume of the plant canopy comprises:
the volume of the ith table was calculated as follows:
the volume of the vertebral body is calculated as follows:
the volume of the plant canopy was calculated as follows:
wherein S isiAnd Si+1The areas of the cross sections corresponding to the lower bottom surface and the upper bottom surface of the ith platform body are respectively, i is the serial number of the platform body and is arranged from the bottom of the plant canopy to the top in sequence, i is 1,2iAnd zi+1Respectively is the z-axis coordinate of the lower bottom surface and the upper bottom surface of the ith platform body; sg-1Is the area of the cross section corresponding to the bottom surface of the vertebral body, zgAnd zg-1Respectively the z-axis coordinates of the base and apex of the vertebral body.
7. A plant canopy volume calculation system based on point cloud data, the system comprising:
the probability density function calculation module is used for projecting the coordinate information of the plant canopy point cloud data to a z axis by taking the three-dimensional coordinate data of the plant canopy point cloud data as an initial value, and calculating a z coordinate probability density function;
the extreme point calculating module is used for calculating a curvature local extreme point of the z coordinate probability density function according to the z coordinate probability density function;
the segmentation module is used for drawing a plane as a cross section for each local extreme point through the point and perpendicular to the z axis; according to a preset number of segmentation layers, segmenting among the local extreme points to obtain a plurality of cross sections perpendicular to the z axis;
the volume calculation module is used for dividing the plant canopy into a cone and a plurality of platforms from top to bottom according to the cross section, and further calculating the volume of the plant canopy;
the method for calculating the volume of the plant canopy comprises the following steps: respectively projecting points in the upper and lower specific distance ranges of each cross section onto the corresponding cross section, and further calculating boundary points of each cross section; connecting the boundary points of each cross section in sequence to obtain a polygonal contour of the corresponding cross section, and further calculating the area of each cross section; dividing the plant canopy into a cone and a plurality of platforms from top to bottom by taking the cross section as an interface; and respectively calculating the volumes of the vertebral body and the platform body according to the area of the cross section, thereby obtaining the volume of the plant canopy.
8. A storage device comprising a memory adapted to store a program, wherein the program is adapted to be loaded and executed by a processor to implement the method of plant canopy volume calculation based on point cloud data of any one of claims 1-6.
9. A processing device, comprising:
a processor adapted to load a program;
a memory adapted to store the program;
characterized in that said program is adapted to be loaded and executed by said processor to implement the plant canopy volume calculation method based on point cloud data of any of claims 1-6.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811571128.8A CN109697731B (en) | 2018-12-21 | 2018-12-21 | Plant canopy volume calculation method and system based on point cloud data and related equipment |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811571128.8A CN109697731B (en) | 2018-12-21 | 2018-12-21 | Plant canopy volume calculation method and system based on point cloud data and related equipment |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109697731A CN109697731A (en) | 2019-04-30 |
| CN109697731B true CN109697731B (en) | 2020-10-27 |
Family
ID=66231869
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811571128.8A Active CN109697731B (en) | 2018-12-21 | 2018-12-21 | Plant canopy volume calculation method and system based on point cloud data and related equipment |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109697731B (en) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110046613B (en) * | 2019-05-16 | 2023-10-24 | 北京农业信息技术研究中心 | Crop canopy in-situ growth phenotype monitoring device and three-dimensional reconstruction method |
| CN110211171A (en) * | 2019-06-18 | 2019-09-06 | 华志微创医疗科技（北京）有限公司 | The method and device of focal area is extracted from medical image |
| CN110889885B (en) * | 2019-10-23 | 2023-05-16 | 中电科新型智慧城市研究院有限公司 | Three-dimensional object volume calculation method based on point cloud slice |
| CN110853092A (en) * | 2019-11-21 | 2020-02-28 | 上海吉七物联网科技有限公司 | Point cloud measurement algorithm based on irregular object |
| CN112837309B (en) * | 2021-03-02 | 2023-10-20 | 华南农业大学 | Fruit tree canopy target recognition device, method, computing equipment and storage medium |
| CN113139569B (en) * | 2021-03-04 | 2022-04-22 | 山东科技大学 | Target classification detection method, device and system |
| CN113610916B (en) * | 2021-06-17 | 2024-04-12 | 同济大学 | Irregular object volume determination method and system based on point cloud data |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104463164A (en) * | 2014-09-03 | 2015-03-25 | 中国科学院遥感与数字地球研究所 | Tree canopy structure information extraction method based on rib method and crown height ratio |
| CN106815847A (en) * | 2017-01-12 | 2017-06-09 | 非凡智慧（宁夏）科技有限公司 | Trees dividing method and single tree extracting method based on laser radar point cloud |
-
2018
- 2018-12-21 CN CN201811571128.8A patent/CN109697731B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104463164A (en) * | 2014-09-03 | 2015-03-25 | 中国科学院遥感与数字地球研究所 | Tree canopy structure information extraction method based on rib method and crown height ratio |
| CN106815847A (en) * | 2017-01-12 | 2017-06-09 | 非凡智慧（宁夏）科技有限公司 | Trees dividing method and single tree extracting method based on laser radar point cloud |
Non-Patent Citations (4)
| Title |
|---|
| Quantification of hidden canopy volume of airborne laser scanning data using a voxel traversal algorithm;Daniel Kükenbrink 等;《Remote Sensing of Environment》;20161022;全文 * |
| 一种基于三维激光扫描***测量树冠体积方法的研究——以油松为例;熊妮娜 等;《北京林业大学学报》;20070831;第62-63页 * |
| 基于地基激光雷达数据的单株玉米三维建模;刘睿 等;《中国农业大学学报》;20140327;全文 * |
| 迭代切片算法在点云曲面拟合中的应用;张崇军 等;《地理空间信息》;20180131;第75页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109697731A (en) | 2019-04-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109697731B (en) | Plant canopy volume calculation method and system based on point cloud data and related equipment | |
| CN108550318B (en) | Map construction method and device | |
| JP4251545B2 (en) | Route planning system for mobile robot | |
| CN107610061B (en) | Edge-preserving point cloud hole repairing method based on two-dimensional projection | |
| CN106023298B (en) | Point cloud Rigid Registration method based on local Poisson curve reestablishing | |
| CN109671155B (en) | Surface grid reconstruction method, system and equipment based on point cloud data | |
| KR101404640B1 (en) | Method and system for image registration | |
| CN111680673B (en) | Method, device and equipment for detecting dynamic object in grid map | |
| CN112347550B (en) | Coupling type indoor three-dimensional semantic graph building and modeling method | |
| CN115564926B (en) | Three-dimensional patch model construction method based on image building structure learning | |
| CN110838115B (en) | Ancient cultural relic three-dimensional model change detection method by contour line extraction and four-dimensional surface fitting | |
| CN110907947B (en) | Real-time loop detection method in mobile robot SLAM problem | |
| JP6840050B2 (en) | Tower tilt detection method and tilt detection device | |
| CN111123242B (en) | Combined calibration method based on laser radar and camera and computer readable storage medium | |
| CN110487274B (en) | SLAM method and system for weak texture scene, navigation vehicle and storage medium | |
| CN110060344B (en) | Prefabricated part overall dimension reverse modeling method based on point cloud data | |
| KR101549155B1 (en) | Method of automatic extraction of building boundary from lidar data | |
| US11734883B2 (en) | Generating mappings of physical spaces from point cloud data | |
| JP6146731B2 (en) | Coordinate correction apparatus, coordinate correction program, and coordinate correction method | |
| CN114396875A (en) | Rectangular parcel volume measurement method based on vertical shooting of depth camera | |
| CN116721228B (en) | Building elevation extraction method and system based on low-density point cloud | |
| KR20090072030A (en) | An implicit geometric regularization of building polygon using lidar data | |
| CN107240133A (en) | Stereoscopic vision mapping model building method | |
| CN116164771A (en) | Mining area road real-time detection method based on semantic map | |
| CN113963114A (en) | Discrete boundary point tracking method based on polygon growth |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |