CN116071530A - Building roof voxelized segmentation method based on airborne laser point cloud - Google Patents

Building roof voxelized segmentation method based on airborne laser point cloud Download PDF

Info

Publication number
CN116071530A
CN116071530A CN202310206932.0A CN202310206932A CN116071530A CN 116071530 A CN116071530 A CN 116071530A CN 202310206932 A CN202310206932 A CN 202310206932A CN 116071530 A CN116071530 A CN 116071530A
Authority
CN
China
Prior art keywords
voxel
coordinate
voxels
points
laser
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.)
Granted
Application number
CN202310206932.0A
Other languages
Chinese (zh)
Other versions
CN116071530B (en
Inventor
张志华
董震
周圣川
杜鹏
张璐琪
门茂林
陈宗强
杨必胜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
QINGDAO INSTITUTE OF SURVEYING AND MAPPING SURVEY
Wuhan University WHU
Original Assignee
QINGDAO INSTITUTE OF SURVEYING AND MAPPING SURVEY
Wuhan University WHU
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by QINGDAO INSTITUTE OF SURVEYING AND MAPPING SURVEY, Wuhan University WHU filed Critical QINGDAO INSTITUTE OF SURVEYING AND MAPPING SURVEY
Priority to CN202310206932.0A priority Critical patent/CN116071530B/en
Publication of CN116071530A publication Critical patent/CN116071530A/en
Application granted granted Critical
Publication of CN116071530B publication Critical patent/CN116071530B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T19/00—Manipulating 3D models or images for computer graphics
    • G06T19/20—Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts

Landscapes

  • Engineering & Computer Science (AREA)
  • Architecture (AREA)
  • Computer Graphics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Generation (AREA)
  • Image Processing (AREA)

Abstract

The invention discloses a building roof voxelization segmentation method based on an airborne laser point cloud. Firstly, building points in an airborne laser point cloud are extracted by using a ground filtering and classifying algorithm, all the building points are segmented into voxels, the coordinate center and normal vector characteristics of the voxels are calculated, then the building points are divided into voxel units, an initial building roof segmentation result is generated by using a region growing algorithm, the voxels are used as nodes, an energy function is built by combining the initial segmentation result and neighborhood voxel characteristic differences, and finally a final building roof segmentation result is obtained by using a graph cutting algorithm. According to the invention, the onboard laser point cloud is utilized to automatically segment the roof structure of the building, the sparse and discrete characteristics of the onboard laser point cloud are overcome, the accurate segmentation of the roof structure of the building in a complex urban scene is realized, the manual intervention can be reduced, and the digital twin, the reconstruction of the urban building model and the like can be served.

Description

Building roof voxelized segmentation method based on airborne laser point cloud
Technical Field
The invention belongs to the technical field of laser scanning data processing, and particularly relates to a building roof voxelization segmentation method based on an airborne laser point cloud.
Background
The three-dimensional laser scanning can scan the earth surface in a top-down mode, can rapidly acquire laser point clouds with geometric coordinates and attributes, and can describe the real world in a digital mode in an active scanning mode, so that the three-dimensional laser scanning becomes an important three-dimensional data acquisition mode and is successfully applied to various geographic information applications. Compared with the two-dimensional image, the three-dimensional laser point cloud can acquire accurate height and geometric structure information of the urban ground object entity. The airborne three-dimensional laser scanning system can rapidly acquire three-dimensional coordinate information of a large-scale city entity, and obtains the city ground object entity structure based on the acquired three-dimensional laser point cloud.
With the advancement of digital twin technology, it is necessary to extract a city building solid model with rich details, and meanwhile, the generation of a refined city building model becomes the foundation of digital substrate construction. In order to better represent urban entities, the live-action three-dimensional model contains the three-dimensional constitutive structure of the entity and the semantic relationships at the component level.
Urban building is complex in composition, and the roofs of high-rise buildings can contain different floor structures; low-rise building roofs may contain roof structures with different slopes. For complex urban building scenes, fine extraction of different structures of the building is crucial for live-action three-dimensional reconstruction. The applications of urban solar energy analysis, disaster emergency, change detection and the like all have requirements on accurate roof structure semantic geometric relationship reconstruction results. Therefore, how to efficiently extract different structures of a building and to effectively divide a building roof in order to construct a finer building model becomes a difficult problem for research.
Disclosure of Invention
Because the point cloud has the characteristic of discreteness, the traditional point cloud clustering or segmentation algorithm is difficult to effectively extract different structures of a building, and a large amount of manual editing is needed. Aiming at the defects of the prior art, the invention provides the building roof voxelized segmentation method based on the airborne laser point cloud, which does not need manual editing, can automatically and accurately segment the building roof, and improves the reliability of automatic modeling.
In order to achieve the above purpose, the technical scheme provided by the invention is a building roof voxelization segmentation method based on an airborne laser point cloud, which comprises the following steps:
step 1, extracting building points in an airborne laser point cloud by using a ground filtering and classifying algorithm;
step 2, dividing all the building points obtained in the step 1 into voxels, and calculating the coordinate center and normal vector of the voxels;
step 3, generating an initial building roof segmentation result by using a region growing algorithm by taking voxels as units;
and 4, constructing an energy function by taking the voxels as nodes and combining the initial segmentation result and the neighborhood voxel characteristic difference, and obtaining a final building roof segmentation result by using a graph segmentation algorithm.
In the step 1, ground points are filtered out by a filtering algorithm for airborne laser point clouds to obtain non-ground points, and then the non-ground points are marked as building points and non-building points point by using a point cloud classification algorithm.
The definition of the airborne laser three-dimensional point cloud is as follows:
Figure SMS_1
（1）
in the method, in the process of the invention,
Figure SMS_2
for the number of laser points in the airborne laser three-dimensional point cloud,/>
Figure SMS_3
is the third point in the three-dimensional point cloudkPoint(s) of (E)>
Figure SMS_4
Is the third point in the three-dimensional point cloudkThe coordinate values of the corresponding space rectangular coordinate system X, Y, Z axes of the points respectively, wherein the space rectangular coordinate system takes any laser point as the origin of coordinates, the X axis points to the east, the Y axis points to the north, and the Z axis forms a right hand system along the normal direction of the ellipsoid of the earth.
In the step 2, the octree structure is used for indexing the building points, the voxels are regularly divided in the point cloud space according to the given voxel size, and the coordinate center and the normal vector of the voxels are calculated by the voxels.
Coordinate center of voxel
Figure SMS_5
The calculation method is as follows:
Figure SMS_6
（2）
in the method, in the process of the invention,
Figure SMS_7
for voxel->
Figure SMS_9
The coordinate centers of all points in the interior are defined by voxel +.>
Figure SMS_12
Vector consisting of spatial coordinate average values of all laser points in the field, < >>
Figure SMS_14
For voxel->
Figure SMS_15
Total number of inner laser spots +.>
Figure SMS_16
For voxel->
Figure SMS_17
Inner firsttCoordinate values of the laser points corresponding to the X axis, < >>
Figure SMS_8
For voxel->
Figure SMS_10
Inner firsttCoordinate values of the laser points corresponding to the Y axis, < >>
Figure SMS_11
For voxel->
Figure SMS_13
Inner firsttCoordinate values corresponding to the laser points on the Z axis.
The normal vector calculation method of the voxels is as follows:
firstly, voxels are utilized
Figure SMS_18
All laser spots contained in the laser beam are used for constructing covariance matrix +.>
Figure SMS_19
The method comprises the following steps:
Figure SMS_20
（3）
in the method, in the process of the invention,
Figure SMS_34
for voxel->
Figure SMS_36
Total number of inner laser spots +.>
Figure SMS_37
For voxel->
Figure SMS_38
Inner firsttVectors composed of spatial coordinates of the individual laser points, +.>
Figure SMS_39
，/>
Figure SMS_40
For voxel->
Figure SMS_41
Inner firsttCoordinate values corresponding to the laser points on the X axis,
Figure SMS_21
for voxel->
Figure SMS_23
Inner firsttCoordinate values of the laser points corresponding to the Y axis, < >>
Figure SMS_25
For voxel->
Figure SMS_28
Inner firsttCoordinate values corresponding to the laser points on the Z axis, < >>
Figure SMS_30
For voxel->
Figure SMS_31
The coordinate centers of all points in the interior are defined by voxel +.>
Figure SMS_33
Vector consisting of spatial coordinate average values of all laser points in the field, < >>
Figure SMS_35
，/>
Figure SMS_22
For voxel->
Figure SMS_24
Coordinate mean value of all laser spots in X-axis, ->
Figure SMS_26
For voxel->
Figure SMS_27
Coordinate mean value of all laser spots in Y-axis, ->
Figure SMS_29
For voxel->
Figure SMS_32
The mean of the coordinates of all laser points in the Z-axis.
Covariance matrix is then formed
Figure SMS_42
Three corresponding eigenvalues are obtained through eigenvalue decomposition
Figure SMS_43
，/>
Figure SMS_44
Eigenvalue->
Figure SMS_45
The corresponding feature vector is voxel->
Figure SMS_46
Corresponding normal vector->
Figure SMS_47
。
And in the step 3, the horizontal and vertical distances between the voxel coordinate center and all the neighborhood voxel coordinate centers and the normal vector included angles between the voxel coordinate center and the neighborhood voxels are calculated respectively, a horizontal distance threshold value, a vertical distance threshold value and a normal vector included angle threshold value of the neighborhood voxel coordinate center are set, the threshold values are used as constraint conditions, the voxels are used as units, and the neighborhood voxels meeting the constraint conditions are grown by using a region growing algorithm, so that an initial building roof segmentation result is obtained.
The horizontal direction distance is calculated by the following formula:
Figure SMS_48
（4）
in the method, in the process of the invention,
Figure SMS_49
for voxel->
Figure SMS_52
Coordinate center and adjacent voxels->
Figure SMS_53
Horizontal distance of coordinate center, +.>
Figure SMS_56
For voxel->
Figure SMS_57
Coordinate mean value of all laser spots in X-axis, ->
Figure SMS_58
For voxel->
Figure SMS_59
Coordinate mean value of all laser spots in Y-axis, ->
Figure SMS_50
For voxel->
Figure SMS_51
Coordinate mean value of all laser spots in X-axis, ->
Figure SMS_54
For voxel->
Figure SMS_55
The average of the coordinates of all laser points in the Y-axis.
The vertical distance is calculated by the following formula:
Figure SMS_60
（5）
in the method, in the process of the invention,
Figure SMS_61
for voxel->
Figure SMS_62
Coordinate center and adjacent voxels->
Figure SMS_63
The vertical distance between the coordinate centers is +.>
Figure SMS_64
Voxel, voxel generation method
Figure SMS_65
Coordinate mean value of all laser spots in Z-axis, < >>
Figure SMS_66
For voxel->
Figure SMS_67
The mean of the coordinates of all laser points in the Z-axis.
Angle of normal vector
Figure SMS_68
The calculation mode of (2) is as follows:
Figure SMS_69
（6）
Figure SMS_70
（7）
in the method, in the process of the invention,
Figure SMS_71
for voxel->
Figure SMS_72
Corresponding normal vector, ++>
Figure SMS_73
For voxel->
Figure SMS_74
Corresponding normal vector, ++>
Figure SMS_75
，
Figure SMS_76
，/>
Figure SMS_77
Is the number of building voxels.
Given a horizontal distance threshold
Figure SMS_79
Vertical distance threshold>
Figure SMS_80
Normal vector angle threshold->
Figure SMS_81
If the condition is satisfied
Figure SMS_82
And->
Figure SMS_83
And->
Figure SMS_84
Then voxel is made +.>
Figure SMS_85
And voxel->
Figure SMS_78
And (5) growing to generate an initial segmentation result.
Moreover, the graph segmentation algorithm in the step 4 is a segmentation method based on graph theory, the core idea is energy optimization, and the graph segmentation algorithm is widely applied to foreground and background segmentation of images in the field of computer vision. Each voxel is subjected to
Figure SMS_86
As a node, optimizing the initial building roof segmentation result by using a graph segmentation algorithm, namely taking the initial segmentation result obtained in the step 3 as an initial label, andconstructing an energy function, taking the minimum energy function as an optimization target, and acquiring a new segmentation label of each voxel by using a graph segmentation algorithm to obtain a final building roof segmentation result; the energy function is calculated as follows, and consists of data items and smooth items: />
Figure SMS_87
（8）
In the method, in the process of the invention,
Figure SMS_89
is a label corresponding to all voxels, < >>
Figure SMS_90
For the number of building voxels +.>
Figure SMS_92
Representing voxel->
Figure SMS_94
The label of (2) is->
Figure SMS_95
The cost value in time, N is the set of the current voxel and all neighborhood voxel group pairs, wherein the neighborhood voxel is voxel +.>
Figure SMS_97
26 neighborhood voxels spatially adjacent in the direction X, Y, Z, respectively, < >>
Figure SMS_99
Index for adjacent voxels, < >>
Figure SMS_88
Is voxel->
Figure SMS_91
And neighborhood voxel->
Figure SMS_93
The labels of (2) are->
Figure SMS_96
And->
Figure SMS_98
Cost value of (a).
Data item
Figure SMS_100
Defined by the initial segmentation result similarity to region growth:
Figure SMS_101
（9）
in the method, in the process of the invention,
Figure SMS_102
for voxel->
Figure SMS_103
Is used for the initial tag of (a),ais voxel->
Figure SMS_104
Label->
Figure SMS_105
With the original tag->
Figure SMS_106
Inconsistent penalty values.
Smoothing terms generated by region growing process
Figure SMS_107
To optimize the over-segmentation phenomenon due to region growth, the smoothing term is therefore defined as:
Figure SMS_108
（10）
Figure SMS_109
（11）
in the method, in the process of the invention,
Figure SMS_111
for voxel->
Figure SMS_112
Coordinate center and neighborhood voxel->
Figure SMS_115
Three-dimensional distance of coordinate center->
Figure SMS_117
For voxel->
Figure SMS_119
And neighborhood voxel->
Figure SMS_120
Is arranged at the angle of the normal vector,bandcpenalty term for tag inconsistency, +.>
Figure SMS_121
For voxel->
Figure SMS_110
Coordinate center and neighborhood voxel->
Figure SMS_113
Three-dimensional distance threshold of coordinate center, +.>
Figure SMS_114
For voxel->
Figure SMS_116
And neighborhood voxel->
Figure SMS_118
A normal vector angle threshold.
And obtaining the optimal labels corresponding to all the voxels through minimizing the energy function, wherein the voxels with the same label finally are the same building roof structure.
Compared with the prior art, the invention has the following advantages:
1) The building roof structure is automatically segmented by utilizing the airborne laser point cloud, the point cloud is divided into voxels and segmented by taking the voxels as units, the sparse and discrete characteristics of the airborne point cloud are overcome, the building roof structure can be accurately segmented, the manual intervention can be reduced, and the efficiency of live-action three-dimensional reconstruction is improved.
2) After dividing the building point cloud into voxels, dividing the building roof structure from thick to thin by using the voxels as units and using a region growing algorithm and a graph cutting algorithm.
Drawings
FIG. 1 is a flow chart of an embodiment of the present invention.
Detailed Description
The invention provides a building roof voxelization segmentation method based on an onboard laser point cloud.
The technical scheme of the invention is further described below with reference to the accompanying drawings and examples.
As shown in fig. 1, the flow of the embodiment of the present invention includes the following steps:
and step 1, extracting building points in the airborne laser point cloud by using a ground filtering and classifying algorithm.
Firstly, filtering ground points of an airborne laser point cloud by using a cloth filtering algorithm to obtain non-ground points, and then marking the non-ground points as building points and non-building points point by using a point cloud classification algorithm.
The definition of the airborne laser three-dimensional point cloud is as follows:
Figure SMS_122
（1）
in the method, in the process of the invention,
Figure SMS_123
for the number of laser points in the airborne laser three-dimensional point cloud,/-for>
Figure SMS_124
Is the third point in the three-dimensional point cloudkPoint(s) of (E)>
Figure SMS_125
Is the third point in the three-dimensional point cloudkThe coordinate values of the corresponding space rectangular coordinate system X, Y, Z axes of the points respectively, wherein the space rectangular coordinate system takes any laser point as the origin of coordinates, the X axis points to the east, the Y axis points to the north, and the Z axis forms a right hand system along the normal direction of the ellipsoid of the earth.
And 2, dividing all the building points obtained in the step 1 into voxels, and calculating the coordinate center and normal vector of the voxels.
And (3) establishing indexes for building points by utilizing the octree structure, dividing voxels in a point cloud space rule according to a given voxel size, and calculating the coordinate center and normal vector of the voxels by voxel.
Coordinate center of voxel
Figure SMS_126
The calculation method is as follows:
Figure SMS_127
（2）
in the method, in the process of the invention,
Figure SMS_129
for voxel->
Figure SMS_131
The coordinate centers of all points in the interior are defined by voxel +.>
Figure SMS_133
Vector consisting of spatial coordinate average values of all laser points in the field, < >>
Figure SMS_135
For voxel->
Figure SMS_136
Total number of inner laser spots +.>
Figure SMS_137
For voxel->
Figure SMS_138
Inner firsttCoordinate values of the laser points corresponding to the X axis, < >>
Figure SMS_128
For voxel->
Figure SMS_130
Inner firsttCoordinate values of the laser points corresponding to the Y axis, < >>
Figure SMS_132
For voxel->
Figure SMS_134
Inner firsttCoordinate values corresponding to the laser points on the Z axis.
The normal vector calculation method of the voxels is as follows:
firstly, voxels are utilized
Figure SMS_139
All laser spots contained in the laser beam are used for constructing covariance matrix +.>
Figure SMS_140
The method comprises the following steps: />
Figure SMS_141
（3）
In the method, in the process of the invention,
Figure SMS_155
for voxel->
Figure SMS_156
Total number of inner laser spots +.>
Figure SMS_158
For voxel->
Figure SMS_159
Inner firsttVectors composed of spatial coordinates of the individual laser points, +.>
Figure SMS_160
，/>
Figure SMS_161
For voxel->
Figure SMS_162
Inner firsttCoordinate values of the laser points corresponding to the X axis, < >>
Figure SMS_142
For voxel->
Figure SMS_145
Inner firsttCoordinate values of the laser points corresponding to the Y axis, < >>
Figure SMS_147
For voxel->
Figure SMS_148
Inner firsttCoordinate values corresponding to the laser points on the Z axis, < >>
Figure SMS_151
For voxel->
Figure SMS_153
The coordinate centers of all points in the interior are defined by voxel +.>
Figure SMS_154
Vector consisting of spatial coordinate average values of all laser points in the field, < >>
Figure SMS_157
，/>
Figure SMS_143
For voxel->
Figure SMS_144
The average of the coordinates of all laser points in the X-axis,
Figure SMS_146
for voxel->
Figure SMS_149
Coordinate mean value of all laser spots in Y-axis, ->
Figure SMS_150
For voxel->
Figure SMS_152
The mean of the coordinates of all laser points in the Z-axis.
Covariance matrix is then formed
Figure SMS_163
Three corresponding eigenvalues are obtained through eigenvalue decomposition
Figure SMS_164
，/>
Figure SMS_165
Eigenvalue->
Figure SMS_166
The corresponding feature vector is voxel->
Figure SMS_167
Corresponding normal vector->
Figure SMS_168
。
And 3, generating an initial building roof segmentation result by using the voxels as units and utilizing a region growing algorithm.
And respectively calculating the horizontal and vertical distances between the voxel coordinate center and all neighbor voxel coordinate centers and the normal vector included angles between the voxel coordinate centers and the neighbor voxels, setting a horizontal distance threshold value, a vertical distance threshold value and a normal vector included angle threshold value of the neighbor voxel coordinate centers, taking the threshold values as constraint conditions, and growing the neighbor voxels meeting the constraint conditions by using a region growing algorithm by taking the voxels as units to obtain an initial building roof segmentation result.
The horizontal direction distance is calculated by the following formula:
Figure SMS_169
（4）
in the method, in the process of the invention,
Figure SMS_171
for voxel->
Figure SMS_172
Coordinate center and adjacent voxels->
Figure SMS_174
Horizontal distance of coordinate center, +.>
Figure SMS_176
For voxel->
Figure SMS_178
Coordinate mean value of all laser spots in X-axis, ->
Figure SMS_179
For voxel->
Figure SMS_180
The average of the coordinates of all laser points in the Y-axis,
Figure SMS_170
for voxel->
Figure SMS_173
Coordinate mean value of all laser spots in X-axis, ->
Figure SMS_175
For voxel->
Figure SMS_177
The average of the coordinates of all laser points in the Y-axis.
The vertical distance is calculated by the following formula:
Figure SMS_181
（5）
in the method, in the process of the invention,
Figure SMS_182
for voxel->
Figure SMS_183
Coordinate center and adjacent voxels->
Figure SMS_184
The vertical distance between the coordinate centers is +.>
Figure SMS_185
Voxel->
Figure SMS_186
Coordinate mean value of all laser spots in Z-axis, < >>
Figure SMS_187
For voxel->
Figure SMS_188
The mean of the coordinates of all laser points in the Z-axis.
Angle of normal vector
Figure SMS_189
The calculation mode of (2) is as follows: />
Figure SMS_190
（6）
Figure SMS_191
（7）
In the method, in the process of the invention,
Figure SMS_192
for voxel->
Figure SMS_193
Corresponding normal vector, ++>
Figure SMS_194
For voxel->
Figure SMS_195
Corresponding normal vector, ++>
Figure SMS_196
，
Figure SMS_197
，/>
Figure SMS_198
Is the number of building voxels.
Given a horizontal distance threshold
Figure SMS_200
Vertical distance threshold>
Figure SMS_201
Normal vector angle threshold->
Figure SMS_204
If the condition is satisfied
Figure SMS_206
And->
Figure SMS_207
And->
Figure SMS_209
Then voxel is made +.>
Figure SMS_210
And voxel->
Figure SMS_199
And (5) growing to generate an initial segmentation result. In this embodiment, the horizontal distance threshold +.>
Figure SMS_202
Is 4m, vertical distance threshold->
Figure SMS_203
Is 4m, normal vector included angle threshold value +.>
Figure SMS_205
Is->
Figure SMS_208
。
And 4, constructing an energy function by taking the voxels as nodes and combining the initial segmentation result and the neighborhood voxel characteristic difference, and obtaining a final building roof segmentation result by using a graph segmentation algorithm.
And taking each voxel as a node, and optimizing the initial building roof segmentation result by using a graph segmentation algorithm. The graph cutting algorithm is a graph theory-based segmentation method, the core idea of which is energy optimization, and the graph cutting algorithm is widely applied to foreground and background segmentation of images in the field of computer vision. And (3) taking the initial segmentation result obtained in the step (3) as an initial label, constructing a data item of an energy function, constructing a smooth item of the energy function by utilizing the geometrical characteristic difference of the neighborhood voxels, taking the minimum energy function as an optimization target, and acquiring a new segmentation label of each voxel by utilizing a graph segmentation algorithm to obtain a final building roof segmentation result.
The graph cut algorithm is treated as a multi-topic problem, i.e., each voxel
Figure SMS_211
All act as one node and are assigned a label in the finite set L. The optimization objective is to give each voxel +.>
Figure SMS_212
Assign a label->
Figure SMS_213
And minimizes the constructed energy function value. The energy function is calculated as follows, and consists of data items and smooth items:
Figure SMS_214
（8）
in the method, in the process of the invention,
Figure SMS_216
is a label corresponding to all voxels, < >>
Figure SMS_217
For the number of building voxels +.>
Figure SMS_220
Representing voxel->
Figure SMS_221
The label of (2) is->
Figure SMS_223
The cost value in time, N is the set of the current voxel and all neighborhood voxel group pairs, wherein the neighborhood voxel is voxel +.>
Figure SMS_225
26 neighborhood voxels spatially adjacent in the direction X, Y, Z, respectively, < >>
Figure SMS_226
Index for adjacent voxels, < >>
Figure SMS_215
Is voxel->
Figure SMS_218
And neighborhood voxel->
Figure SMS_219
The labels of (2) are->
Figure SMS_222
And->
Figure SMS_224
Cost value of (a).
The data item is defined by an initial segmentation result similarity to the region growth:
Figure SMS_227
（9）/>
in the method, in the process of the invention,
Figure SMS_228
for voxel->
Figure SMS_229
Is used for the initial tag of (a),ais voxel->
Figure SMS_230
Label->
Figure SMS_231
With the original tag->
Figure SMS_232
Inconsistent penalty values, this embodimentaThe value was taken as 5.
The smoothing term generated in the region growing process is to optimize the over-segmentation phenomenon caused by the region growing, so the smoothing term is defined as:
Figure SMS_233
（10）
Figure SMS_234
（11）
in the method, in the process of the invention,
Figure SMS_236
for voxel->
Figure SMS_238
Coordinate center and neighborhood voxel->
Figure SMS_240
Three-dimensional distance of coordinate center->
Figure SMS_242
Is a voxel
Figure SMS_244
And neighborhood voxel->
Figure SMS_246
Is arranged at the angle of the normal vector,bandcpenalty term for tag inconsistency, +.>
Figure SMS_248
For voxel->
Figure SMS_235
Coordinate center and neighborhood voxel->
Figure SMS_237
Three-dimensional distance threshold of coordinate center, +.>
Figure SMS_239
For voxel->
Figure SMS_241
And neighborhood voxel->
Figure SMS_243
Normal vector angle threshold, in this embodiment, three-dimensional distance threshold +.>
Figure SMS_245
6m, normal vector angle threshold +.>
Figure SMS_247
Is->
Figure SMS_249
Penalty valuebAndcand are all 7.
By using
Figure SMS_250
The expansion algorithm minimizes the energy function to obtain the best label corresponding to all voxels, and finally the voxels with the same label are the same building roof structure.
The specific embodiments described herein are offered by way of example only to illustrate the spirit of the invention. Those skilled in the art may make various modifications or additions to the described embodiments or substitutions thereof without departing from the spirit of the invention or exceeding the scope of the invention as defined in the accompanying claims.

Claims (9)

1. The building roof voxelization segmentation method based on the airborne laser point cloud is characterized by comprising the following steps of:
step 1, extracting building points in an airborne laser point cloud by using a ground filtering and classifying algorithm;
step 2, dividing all the building points obtained in the step 1 into voxels, and calculating the coordinate center and normal vector of the voxels;
step 3, generating an initial building roof segmentation result by using a region growing algorithm by taking voxels as units;
step 4, using voxels as nodes, constructing an energy function by combining the initial segmentation result and the neighborhood voxel characteristic difference, and obtaining a final building roof segmentation result by using a graph segmentation algorithm;
the graph cutting algorithm is a graph theory-based segmentation method, the core idea of which is energy optimization, and the graph cutting algorithm is widely applied to foreground and background segmentation of images in the field of computer vision.
2. The method for voxel segmentation of building roof based on airborne laser point cloud as set forth in claim 1, wherein: in the step 1, firstly, ground points are filtered out by using a filtering algorithm for airborne laser point clouds to obtain non-ground points, then, the non-ground points are marked as building points and non-building points point by using a point cloud classification algorithm, and the definition of the airborne laser three-dimensional point clouds is as follows:
Figure QLYQS_1
（1）
in the method, in the process of the invention,
Figure QLYQS_2
for the number of laser points in the airborne laser three-dimensional point cloud,/-for>
Figure QLYQS_3
Is the third point in the three-dimensional point cloudkA point of the light-emitting diode is located,
Figure QLYQS_4
is the third point in the three-dimensional point cloudkThe coordinate values of the corresponding space rectangular coordinate system X, Y, Z axes of the points respectively, wherein the space rectangular coordinate system takes any laser point as the origin of coordinates, the X axis points to the east, the Y axis points to the north, and the Z axis forms a right hand system along the normal direction of the ellipsoid of the earth.
3. The method for voxel segmentation of building roof based on airborne laser point cloud as set forth in claim 1, wherein: in step 2, building points are indexed by utilizing an octree structure, voxels are divided in a point cloud space rule according to given voxel sizes, coordinate centers and normal vectors of the voxels are calculated from voxel to voxel, and the coordinate centers of the voxels are calculated
Figure QLYQS_5
The calculation method is as follows:
Figure QLYQS_6
（2）
in the method, in the process of the invention,
Figure QLYQS_8
for voxel->
Figure QLYQS_10
The coordinate centers of all points in the interior are defined by voxel +.>
Figure QLYQS_11
Vector consisting of spatial coordinate average values of all laser points in the field, < >>
Figure QLYQS_14
Is a voxel/>
Figure QLYQS_15
Total number of inner laser spots +.>
Figure QLYQS_16
For voxel->
Figure QLYQS_17
Inner firsttCoordinate values of the laser points corresponding to the X axis, < >>
Figure QLYQS_7
For voxel->
Figure QLYQS_9
Inner firsttCoordinate values of the laser points corresponding to the Y axis, < >>
Figure QLYQS_12
For voxel->
Figure QLYQS_13
Inner firsttCoordinate values corresponding to the laser points on the Z axis.
4. A method for voxelized segmentation of building roofs based on an on-board laser point cloud as set forth in claim 3, wherein: the normal vector calculation method of the voxels in the step 2 is as follows: firstly, voxels are utilized
Figure QLYQS_18
All laser spots contained in the laser beam are used for constructing covariance matrix +.>
Figure QLYQS_19
The method comprises the following steps:
Figure QLYQS_20
（3）
in the method, in the process of the invention,
Figure QLYQS_32
for voxel->
Figure QLYQS_35
Total number of inner laser spots +.>
Figure QLYQS_37
For voxel->
Figure QLYQS_38
Inner firsttVectors composed of spatial coordinates of the individual laser points, +.>
Figure QLYQS_39
，/>
Figure QLYQS_40
For voxel->
Figure QLYQS_41
Inner firsttCoordinate values of the laser points corresponding to the X axis, < >>
Figure QLYQS_21
For voxel->
Figure QLYQS_24
Inner firsttCoordinate values of the laser points corresponding to the Y axis, < >>
Figure QLYQS_26
For voxel->
Figure QLYQS_28
Inner firsttCoordinate values corresponding to the laser points on the Z axis, < >>
Figure QLYQS_30
For voxel->
Figure QLYQS_33
All points inThe coordinate center is composed of voxels->
Figure QLYQS_34
Vector consisting of spatial coordinate average values of all laser points in the field, < >>
Figure QLYQS_36
，/>
Figure QLYQS_22
For voxel->
Figure QLYQS_23
Coordinate mean value of all laser spots in X-axis, ->
Figure QLYQS_25
For voxel->
Figure QLYQS_27
Coordinate mean value of all laser spots in Y-axis, ->
Figure QLYQS_29
For voxel->
Figure QLYQS_31
The average value of the coordinates of all laser points in the Z axis;
covariance matrix is then formed
Figure QLYQS_42
The corresponding three eigenvalues +.>
Figure QLYQS_43
，
Figure QLYQS_44
Eigenvalue->
Figure QLYQS_45
The corresponding feature vector is voxel->
Figure QLYQS_46
Corresponding normal vector->
Figure QLYQS_47
。
5. The method for voxel segmentation of building roof based on airborne laser point cloud as set forth in claim 4, wherein: in the step 3, respectively calculating the horizontal and vertical distances between the voxel coordinate center and all neighborhood voxel coordinate centers and the normal vector included angles between the voxel coordinate centers and the neighborhood voxels, setting a horizontal distance threshold value, a vertical distance threshold value and a normal vector included angle threshold value of the neighborhood voxel coordinate centers, taking the threshold values as constraint conditions, and growing the neighborhood voxels meeting the constraint conditions by using a region growing algorithm by taking the voxels as units to obtain an initial building roof segmentation result; the horizontal direction distance is calculated by the following formula:
Figure QLYQS_48
（4）
in the method, in the process of the invention,
Figure QLYQS_50
for voxel->
Figure QLYQS_51
Coordinate center and adjacent voxels->
Figure QLYQS_53
Horizontal distance of coordinate center, +.>
Figure QLYQS_55
Is a voxel
Figure QLYQS_57
Coordinate mean value of all laser spots in X-axis, ->
Figure QLYQS_58
For voxel->
Figure QLYQS_59
Coordinate mean value of all laser spots in Y-axis, ->
Figure QLYQS_49
For voxel->
Figure QLYQS_52
Coordinate mean value of all laser spots in X-axis, ->
Figure QLYQS_54
For voxel->
Figure QLYQS_56
The average value of the coordinates of all laser points in the Y axis;
the vertical distance is calculated by the following formula:
Figure QLYQS_60
（5）
in the method, in the process of the invention,
Figure QLYQS_61
for voxel->
Figure QLYQS_62
Coordinate center and adjacent voxels->
Figure QLYQS_63
The vertical distance between the coordinate centers is +.>
Figure QLYQS_64
Voxel->
Figure QLYQS_65
Coordinate mean value of all laser spots in Z-axis, < >>
Figure QLYQS_66
For voxel->
Figure QLYQS_67
The mean of the coordinates of all laser points in the Z-axis.
6. The method for voxel segmentation of building roof based on airborne laser point cloud as set forth in claim 5, wherein: normal vector included angle in step 3
Figure QLYQS_68
The calculation mode of (2) is as follows:
Figure QLYQS_69
（6）
Figure QLYQS_70
（7）
in the method, in the process of the invention,
Figure QLYQS_71
for voxel->
Figure QLYQS_72
Corresponding normal vector, ++>
Figure QLYQS_73
For voxel->
Figure QLYQS_74
Corresponding normal vector, ++>
Figure QLYQS_75
，/>
Figure QLYQS_76
，/>
Figure QLYQS_77
Number of voxels for the building;
given a horizontal distance threshold
Figure QLYQS_79
Vertical distance threshold>
Figure QLYQS_80
Normal vector angle threshold->
Figure QLYQS_81
If the condition is satisfied
Figure QLYQS_82
And->
Figure QLYQS_83
And->
Figure QLYQS_84
Then voxel is made +.>
Figure QLYQS_85
And voxel->
Figure QLYQS_78
And (5) growing to generate an initial segmentation result.
7. The method for voxel segmentation of building roof based on airborne laser point cloud as set forth in claim 6, wherein: step 4, each voxel is processed
Figure QLYQS_86
As a node, optimizing the roof segmentation result of the initial building by using a graph segmentation algorithm, namely taking the initial segmentation result obtained in the step 3 as an initial label, constructing an energy function, taking the minimum energy function as an optimization target, and obtaining a new segmentation label of each voxel by using the graph segmentation algorithm to obtain a final buildingRoof segmentation results; the energy function is calculated as follows, and consists of data items and smooth items:
Figure QLYQS_87
（8）
in the method, in the process of the invention,
Figure QLYQS_89
is a label corresponding to all voxels, < >>
Figure QLYQS_91
For the number of building voxels +.>
Figure QLYQS_92
Representing voxel->
Figure QLYQS_95
The label of (2) is->
Figure QLYQS_97
The cost value in time, N is the set of the current voxel and all neighborhood voxel group pairs, wherein the neighborhood voxel is voxel +.>
Figure QLYQS_98
26 neighborhood voxels spatially adjacent in the direction X, Y, Z, respectively, < >>
Figure QLYQS_99
Index for adjacent voxels, < >>
Figure QLYQS_88
Is voxel->
Figure QLYQS_90
And neighborhood voxel->
Figure QLYQS_93
The labels of (2) are->
Figure QLYQS_94
And->
Figure QLYQS_96
Cost value of (a).
8. The method for voxel segmentation of building roof based on airborne laser point cloud as set forth in claim 7, wherein: data item in step 4
Figure QLYQS_100
Defined by the initial segmentation result similarity to region growth:
Figure QLYQS_101
（9）
in the method, in the process of the invention,
Figure QLYQS_102
for voxel->
Figure QLYQS_103
Is used for the initial tag of (a),ais voxel->
Figure QLYQS_104
Label->
Figure QLYQS_105
With the original tag->
Figure QLYQS_106
Inconsistent penalty values.
9. The method for voxelized segmentation of building roofs based on onboard laser point clouds as set forth in claim 8, wherein: smoothing term generated in step 4 in the region growing process
Figure QLYQS_107
Is defined as：
Figure QLYQS_108
（10）
Figure QLYQS_109
（11）
In the method, in the process of the invention,
Figure QLYQS_111
for voxel->
Figure QLYQS_112
Coordinate center and neighborhood voxel->
Figure QLYQS_115
Three-dimensional distance of coordinate center->
Figure QLYQS_116
For voxel->
Figure QLYQS_118
And neighborhood voxel->
Figure QLYQS_120
Is arranged at the angle of the normal vector,bandcpenalty term for tag inconsistency, +.>
Figure QLYQS_121
For voxel->
Figure QLYQS_110
Coordinate center and neighborhood voxel->
Figure QLYQS_113
Three-dimensional distance threshold of coordinate center, +.>
Figure QLYQS_114
For voxel->
Figure QLYQS_117
And neighborhood voxel->
Figure QLYQS_119
A normal vector angle threshold;
and obtaining the optimal labels corresponding to all the voxels through minimizing the energy function, wherein the voxels with the same label finally are the same building roof structure.
CN202310206932.0A 2023-03-07 2023-03-07 Building roof voxelized segmentation method based on airborne laser point cloud Active CN116071530B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310206932.0A CN116071530B (en) 2023-03-07 2023-03-07 Building roof voxelized segmentation method based on airborne laser point cloud

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310206932.0A CN116071530B (en) 2023-03-07 2023-03-07 Building roof voxelized segmentation method based on airborne laser point cloud

Publications (2)

Publication Number Publication Date
CN116071530A true CN116071530A (en) 2023-05-05
CN116071530B CN116071530B (en) 2023-07-28

Family

ID=86180314

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310206932.0A Active CN116071530B (en) 2023-03-07 2023-03-07 Building roof voxelized segmentation method based on airborne laser point cloud

Country Status (1)

Country Link
CN (1) CN116071530B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117576144A (en) * 2024-01-15 2024-02-20 湖北工业大学 Laser point cloud power line extraction method and device and electronic equipment

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090310867A1 (en) * 2008-06-12 2009-12-17 Bogdan Calin Mihai Matei Building segmentation for densely built urban regions using aerial lidar data
US20100171740A1 (en) * 2008-03-28 2010-07-08 Schlumberger Technology Corporation Visualizing region growing in three dimensional voxel volumes
US20130057653A1 (en) * 2011-09-06 2013-03-07 Electronics And Telecommunications Research Institute Apparatus and method for rendering point cloud using voxel grid
CN109242855A (en) * 2018-07-19 2019-01-18 中国科学院自动化研究所 Roof dividing method, system and equipment based on Three-dimensional Multi-resolution statistical information
CN109685821A (en) * 2018-12-26 2019-04-26 中国科学院大学 Region growing 3D rock mass point cloud plane extracting method based on high quality voxel
WO2021097618A1 (en) * 2019-11-18 2021-05-27 深圳市大疆创新科技有限公司 Point cloud segmentation method and system, and computer storage medium
CN113205529A (en) * 2021-04-19 2021-08-03 武汉大学 Method for segmenting top surface of building based on airborne LiDAR point cloud
CN114049395A (en) * 2021-03-31 2022-02-15 武汉大学 Airborne point cloud assisted vehicle-mounted laser point cloud position precision improving method
CN114463338A (en) * 2022-01-07 2022-05-10 武汉大学 Automatic building laser foot point extraction method based on graph cutting and post-processing
CN114758323A (en) * 2022-06-13 2022-07-15 青岛市勘察测绘研究院 Urban road sign extraction method based on vehicle-mounted laser point cloud
CN114764871A (en) * 2022-06-15 2022-07-19 青岛市勘察测绘研究院 Urban building attribute extraction method based on airborne laser point cloud

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100171740A1 (en) * 2008-03-28 2010-07-08 Schlumberger Technology Corporation Visualizing region growing in three dimensional voxel volumes
US20090310867A1 (en) * 2008-06-12 2009-12-17 Bogdan Calin Mihai Matei Building segmentation for densely built urban regions using aerial lidar data
US20130057653A1 (en) * 2011-09-06 2013-03-07 Electronics And Telecommunications Research Institute Apparatus and method for rendering point cloud using voxel grid
CN109242855A (en) * 2018-07-19 2019-01-18 中国科学院自动化研究所 Roof dividing method, system and equipment based on Three-dimensional Multi-resolution statistical information
CN109685821A (en) * 2018-12-26 2019-04-26 中国科学院大学 Region growing 3D rock mass point cloud plane extracting method based on high quality voxel
WO2021097618A1 (en) * 2019-11-18 2021-05-27 深圳市大疆创新科技有限公司 Point cloud segmentation method and system, and computer storage medium
CN114049395A (en) * 2021-03-31 2022-02-15 武汉大学 Airborne point cloud assisted vehicle-mounted laser point cloud position precision improving method
CN113205529A (en) * 2021-04-19 2021-08-03 武汉大学 Method for segmenting top surface of building based on airborne LiDAR point cloud
CN114463338A (en) * 2022-01-07 2022-05-10 武汉大学 Automatic building laser foot point extraction method based on graph cutting and post-processing
CN114758323A (en) * 2022-06-13 2022-07-15 青岛市勘察测绘研究院 Urban road sign extraction method based on vehicle-mounted laser point cloud
CN114764871A (en) * 2022-06-15 2022-07-19 青岛市勘察测绘研究院 Urban building attribute extraction method based on airborne laser point cloud

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
YUSHENG XU 等: "Segmentation of building roofs from airborne LiDAR point clouds using robust voxel-based region growing", 《REMOTE SENSING LETTERS》, pages 1061 - 1071 *
郑特 等: "基于图割算法的摄影测量点云面向对象分类方法", 《测绘工程》, pages 16 - 19 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117576144A (en) * 2024-01-15 2024-02-20 湖北工业大学 Laser point cloud power line extraction method and device and electronic equipment
CN117576144B (en) * 2024-01-15 2024-03-29 湖北工业大学 Laser point cloud power line extraction method and device and electronic equipment

Also Published As

Publication number Publication date
CN116071530B (en) 2023-07-28

Similar Documents

Publication Publication Date Title
CN110570428B (en) Method and system for dividing building roof sheet from large-scale image dense matching point cloud
CN112070769B (en) Layered point cloud segmentation method based on DBSCAN
CN111915730B (en) Method and system for automatically generating indoor three-dimensional model by taking semantic slave point cloud into consideration
CN109961440B (en) Three-dimensional laser radar point cloud target segmentation method based on depth map
CN107292276B (en) Vehicle-mounted point cloud clustering method and system
CN106780524B (en) Automatic extraction method for three-dimensional point cloud road boundary
Sun et al. Aerial 3D building detection and modeling from airborne LiDAR point clouds
Vosselman et al. Recognising structure in laser scanner point clouds
CN111815776A (en) Three-dimensional building fine geometric reconstruction method integrating airborne and vehicle-mounted three-dimensional laser point clouds and streetscape images
CN111598916A (en) Preparation method of indoor occupancy grid map based on RGB-D information
CN112288857A (en) Robot semantic map object recognition method based on deep learning
CN112164145B (en) Method for rapidly extracting indoor three-dimensional line segment structure based on point cloud data
CN114119902A (en) Building extraction method based on unmanned aerial vehicle inclined three-dimensional model
CN117132915B (en) Power transmission line tree obstacle hidden danger analysis method based on automatic classification of point cloud
CN116071530B (en) Building roof voxelized segmentation method based on airborne laser point cloud
CN116258857A (en) Outdoor tree-oriented laser point cloud segmentation and extraction method
CN114463396B (en) Point cloud registration method utilizing plane shape and topological graph voting
Guinard et al. Piecewise-planar approximation of large 3d data as graph-structured optimization
CN114463338A (en) Automatic building laser foot point extraction method based on graph cutting and post-processing
Sun et al. Automated segmentation of LiDAR point clouds for building rooftop extraction
CN113985435A (en) Mapping method and system fusing multiple laser radars
CN113409332A (en) Building plane segmentation method based on three-dimensional point cloud
CN113536959A (en) Dynamic obstacle detection method based on stereoscopic vision
Rothermel et al. Fast and robust generation of semantic urban terrain models from UAV video streams
Lei et al. Automatic identification of street trees with improved RandLA-Net and accurate calculation of shading area with density-based iterative α-shape

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