CN110570428B

CN110570428B - Method and system for dividing building roof sheet from large-scale image dense matching point cloud

Info

Publication number: CN110570428B
Application number: CN201910734783.9A
Authority: CN
Inventors: 王瑞胜; 彭飞宇; 朱正荣; 蒲冰鑫; 张新梅; 钟若飞
Original assignee: Zhejiang Hexin Geographic Information Technology Co ltd
Current assignee: Zhejiang Hexin Geographic Information Technology Co ltd
Priority date: 2019-08-09
Filing date: 2019-08-09
Publication date: 2023-07-07
Anticipated expiration: 2039-08-09
Also published as: CN110570428A

Abstract

The invention discloses a method and a system for dividing a building roof panel from a large-scale image densely-matched point cloud, which aim to solve the problem of how to accurately obtain the building roof three-dimensional panel from a color three-dimensional point cloud generated by an inclined image. A method of segmenting building roof sheeting from a large scale image dense matching point cloud, comprising the steps of: ground filtration and vegetation filtration; removing the ground and vegetation in the image point cloud; building singulation; eliminating non-building points and realizing building clustering; removing the building elevation; cutting, repairing and refining the roof sheet.

Description

Method and system for dividing building roof sheet from large-scale image dense matching point cloud

Technical Field

The invention relates to the technical field of image application, in particular to a method and a system for dividing a building roof sheet from a large-scale image dense matching point cloud.

Background

Cities are changing day-to-day as an important body of spatiotemporal information. The three-dimensional model of the urban building is used as the visual expression of information in a series of urban spaces such as urban topography, buildings on the ground and the like, displays various information such as the geometric shape, the attribute, the spatial position, the texture and the like of the target, and is the basis of the three-dimensional geographic information system of the city. The rapid extraction and detection of the building change have important roles in the aspects of GIS database updating, land utilization, digital city and the like.

The method takes the three-dimensional color point cloud acquired by the oblique photogrammetry system as a data source to study the roof sheet segmentation technology of the building. The following describes in detail the tilt photogrammetry system, the on-board LiDAR based building rooftop extraction technique, respectively.

The oblique photogrammetry system is a new technology developed on the basis of the traditional photogrammetry technology by combining with the computer vision, and can acquire ground images of multiple angles at the same time, and acquire rich textures of the top image and the side surface of a building. The technology has the advantages that the real situation and the form of the ground object can be comprehensively and comprehensively presented, the high-precision side texture information of the ground object, especially the building, can be captured, the ground object in the image can be positioned and modeled, and a real city three-dimensional model can be constructed, so that the application field of the technology becomes wider and wider.

Currently, the tilt imaging system is basically composed of a "1+n" camera, i.e. a down-view camera and N tilt cameras, and from the perspective of three-dimensional modeling in city, the "1+n" tilt imaging system can acquire images with complete coverage. Currently, the most common is 4 tilt cameras, namely, 2 front-back, left-right. The tilt imaging system basically integrates a high-precision POS system, and records the position and posture data of the camera at the moment of image exposure while image acquisition is carried out.

In order to further improve the completeness of image coverage of the oblique imaging system, the adopted technical strategy is to increase the number of oblique cameras, such as OCTOBLIQUE MIDAS of a Track Air formula, to 8, and shooting by a down-view camera and eight oblique cameras with an angle of inclination of 45 degrees in a 360-degree panorama, wherein the additionally added four cameras can create favorable conditions for covering the shot area of the image without dead angles in all directions.

In summary, the multi-view tilt camera is generally mounted on an aircraft, and acquires ground object images at multiple angles. From the technical characteristics of multi-view oblique photography, the following 4 characteristics are mainly provided:

(1) The images shot by the five lenses at the same moment are not overlapped;

(2) In the whole area, the overlapping relation between the observation image and the vertical image is complex and unknown;

(3) A rotation of about 180 degrees exists between the vertical image and the rear view image of the same course having an overlapping region;

(4) In the absence of intersecting routes, there is a rotation of about + -90 DEG between the vertical image and the left-right view image with the overlapping area.

The image-based dense point cloud three-dimensional reconstruction technology, also called multi-view stereo (MVS), uses a scaled image obtained from an SFM algorithm, that is, parameters of a camera corresponding to the image have been obtained through calculation or other methods. Currently, post-processing commercial software mainly includes ContextCapture software produced by Acute3D company in France, altizure platform established by Diluon of university of hong Kong science and technology, and the like. The dense point cloud three-dimensional reconstruction technology based on the image mainly comprises the following steps:

(1) Detecting and matching the feature points;

(2) Solving camera parameters; usually, a camera solving method in the SFM sense is used, under the condition of inaccurate initialization, accurate camera parameters are obtained by utilizing bundle adjustment constraint, and a sparse point cloud model is obtained;

(3) Three-dimensional reconstruction of dense point clouds; at present, a PMVS algorithm is mainly popular, namely, a sparse point cloud and camera parameters obtained by SFM are utilized, and a three-dimensional point cloud model under dense matching is generated through a three-dimensional reconstruction algorithm generated by PMVS based on a patch.

(4) Mesh hole patching and adding texture mapping.

The oblique photogrammetry technology combined with computer vision breaks through the limitations of remote sensing technology, traditional photogrammetry and mapping industry, and has wide development prospect. The technology has great potential in market application in the fields of photogrammetry and remote sensing in China, and has great promotion effect on the construction of intelligent city and digital city. In the method of the invention, oblique photography is the basis of research, and the dense colored three-dimensional point cloud obtained by the method is the object of building roof contour extraction.

In general, building roof extraction technology based on airborne LiDAR is mainly 2 steps: building singulation and roof 3D patch generation. In building singulation, two concepts are commonly used. 1) Identifying the building directly from the original point cloud; 2) And filtering out the ground and vegetation respectively to obtain building point clouds. And the roof sheet extraction technology is carried out according to the removal of the building vertical face and the sheet segmentation.

Disclosure of Invention

In order to solve the technical problems, the invention aims to provide a method and a system for dividing building roof tiles from large-scale image densely matched point clouds, which are used for solving the problem of accurately acquiring the building roof three-dimensional tiles from color three-dimensional point clouds generated by oblique images.

In order to achieve the above purpose, the present invention adopts the following technical scheme:

a method of segmenting building roof sheeting from a large scale image dense matching point cloud, comprising the steps of:

1) Ground filtration and vegetation filtration; removing the ground and vegetation in the image point cloud;

2) Building singulation; eliminating non-building points and realizing building clustering;

3) Removing the building elevation;

4) Cutting, repairing and refining the roof sheet.

Preferably, the ground filtering is performed, and the ground points in the point cloud data are separated from the non-ground points by using a cloth filtering algorithm, so as to obtain the non-ground points.

Preferably, the specific implementation process of the ground filtration is as follows:

1) Turning the original laser point cloud with outliers removed by 180 degrees;

2) Initializing a distribution grid; determining the total number of cloth particles by using the defined grid resolution; the initial position of the cloth is placed above the highest point of the point cloud;

3) Projecting the laser point cloud and the grid particles to a horizontal plane, finding the nearest neighbor laser point of each particle, and taking the height Cheng Ji of the laser points as an intersection point elevation value;

4) For each particle, calculating its position due to gravity; if its elevation value is greater than the intersection elevation value, it can continue to move; if the height value of the particle is smaller than or equal to the height value of the intersection point, the height value of the particle is set as the height value of the intersection point and is marked as immovable;

5) For each grid particle, calculating a positional shift due to an interaction force between the particles;

6) Iteratively performing steps 4) to 5) above until the maximum elevation variance of all particles is sufficiently small or exceeds a defined maximum number of iterations;

7) Calculating the distance between the particles and the point cloud;

8) Separating ground and non-ground points; for any one laser point, if the distance from the corresponding grid particle point is less than a threshold value, identifying the laser point as a ground point; otherwise, it is a non-ground point.

Preferably, the specific implementation process of the ground filtration further comprises a post-processing step: searching for the four neighboring particles of each movable particle, and if an immovable particle is found, comparing the two to the elevation of the corresponding spot; if the minimum elevation difference is less than the threshold value, the movable particles continue to move until set to be immovable; this step is iteratively performed until all movable particles are traversed.

Preferably, the vegetation filtration utilizes a green to blue ratio G _r And identifying green or near-green points by using the vegetation coefficient index VI, optimizing the identification point set by using a variation detection algorithm based on a normal vector, and removing vegetation points.

Preferably, the green signal ratio is a ratio of the DN value of the green band to the sum of the DN values of RGB; the green-to-signal ratio and vegetation coefficient are defined as follows:

wherein: VI and G _r Respectively vegetation coefficient and green letter ratio; r, G, B is the 3 color channel values of the point cloud, R ', G', B 'are the proportions of R, G, B, and R', G ', B' are defined as follows:

preferably, the building is subjected to the building singulation, the connected domain analysis is performed by using an European clustering algorithm to obtain a preliminary building singulation result, and then the constraint term based on the elevation threshold is used for further optimization.

Preferably, removing the building elevation, and removing the building elevation by using the point characteristic f based on the normal vector to obtain a building roof point cloud; the definition of the point feature f is shown in the following formula:

f＝1-|n _p ·e _z |

wherein n is _p Is the normal vector at point p, e _z Is a unit column vector (0, 1).

Preferably, the division of the roof sheet refers to the division of the roof sheet by utilizing a region growing algorithm based on the point spacing for the roof point cloud of a single building; the repairing of the roof sheet refers to the steps of repairing holes by using a sheet repairing algorithm based on PCA on output of the segmented roof sheet and flattening the sheet; the refinement of the roof sheet refers to the output of the divided roof sheet, and the duplication removal of the roof sheet is realized by utilizing a two-dimensional grid.

A system for partitioning building roof sheeting from a large-scale image dense matching point cloud, using a method as described above, comprising:

the first module is used for separating the ground points by using a cloth filtering algorithm to obtain non-ground points;

the second module is used for obtaining initial vegetation points by utilizing the green signal ratio and the vegetation coefficient, and detecting further constraint vegetation points based on the change of the normal vector;

the third module is used for executing European clustering algorithm and elevation threshold constraint so as to realize building singulation;

a fourth module for realizing a building elevation elimination algorithm and roof sheet segmentation to obtain an initial roof sheet;

and a fifth module, configured to implement a patch repair algorithm based on PCA and patch de-duplication based on two-dimensional grid, to obtain a flat roof patch.

By means of some ideas of LiDAR point cloud processing, the invention provides an algorithm framework for fully automatically extracting three-dimensional contour lines of a building from an image point cloud. Because the image point cloud is much noisier, the algorithm is more complex than the LiDAR point cloud. According to the flow, the framework mainly comprises three parts: building singulation, roof sheet segmentation and boundary line generation. Wherein each part consists of 2-3 basic algorithms. The contents of each section are described in detail below.

1. Building singulation

Building singulation is an essential step in order to achieve accurate extraction of building three-dimensional contours from a large-scale image point cloud. Because the roof contour line is extracted by taking a single building as a unit, the contour line is extracted on the basis. By means of the idea that LiDAR point clouds separate buildings from bottom to top, the method realizes building singulation by separating the ground and vegetation in sequence.

1.1 floor filtration

The image point cloud is similar to the on-board laser point cloud in terms of data quality analysis, so that ground removal can be performed using ground identification algorithms in the on-board laser point cloud. The cloth filtering algorithm proposed by Wuming Zhang in 2016 is used herein, and the paper is detailed in Wuming Z, jianbo Q, peng W, et al, an Easy-to-Use Airborne LiDAR Data Filtering Method Based on Cloth Simulation [ J ], remote Sensing,2016,8 (6): 501-. The algorithm has less parameter setting, strong algorithm generalization capability and stable ground filtering effect.

The cloth filtering algorithm simulates a simple physical process. It is assumed that a piece of cloth falls off in the air due to gravity. If the cloth is soft enough to adhere to the object surface, the final shape of the cloth is a digital surface model of the scene. The point cloud scene is now flipped 180 degrees with the ground portion in the point cloud above. In contrast, where the cloth has a certain stiffness, the final shape of the cloth is a digital terrain model, i.e. the ground. The final shape of the cloth can be determined by analyzing the interaction of the particles of the cloth and the corresponding laser point cloud, thereby separating the ground points from the non-ground points. Based on this idea, the algorithm uses a cloth simulation technique in three-dimensional computer graphics, i.e. by computer programming to simulate a cloth. The cloth here is a grid formed by a large number of particles of fixed mass but without space size connected. The vertices of the mesh are these physical particles and the sides of the mesh are a "spring" that obeys the holk's law.

The specific implementation process of the algorithm is as follows:

1) Turning the original laser point cloud with outliers removed by 180 degrees;

2) Initializing the distribution grid. Determining the total number of cloth particles according to the user-defined grid resolution; the initial position of the cloth is placed above the highest point of the point cloud;

3) The laser point cloud and grid particles are projected to a horizontal plane. Finding the nearest neighbor laser point of each particle, and taking the height Cheng Ji of the laser point as an intersection point elevation value;

4) For each particle, its position after the effect of gravity is calculated. If its elevation value is greater than the intersection elevation value, it can continue to move; if the height value of the particle is smaller than or equal to the height value of the intersection point, the height value of the particle is set as the height value of the intersection point and is marked as immovable;

6) Iteratively performing 4) -5) until the maximum elevation variance of all particles is sufficiently small or exceeds a user-defined maximum number of iterations;

7) Calculating the distance between the particles and the point cloud;

8) Separating ground and non-ground points. For any one laser point, if the distance from the corresponding grid particle point is less than a threshold value, identifying the laser point as a ground point; otherwise, the non-ground point is obtained;

In order to adapt the algorithm to terrain with large gradient changes, a post-processing algorithm is added. The principle is as follows: the four neighboring particles of each movable particle are searched, and if an immovable particle is found, the heights of the two to the corresponding spot are compared. If the minimum elevation difference is less than the threshold value, the movable particles continue to move until set to be immovable. This step is iteratively performed until all movable particles are traversed.

1.2 Vegetation Filter

The present algorithm does not employ the idea of a conventional laser point cloud. From the data content analysis, the data source is a colored three-dimensional point cloud with correct RGB information; since vegetation is typically green or near-green, the green band DN (Digital Number) value of vegetation is typically higher than the red and blue band DN values. In combination with the two points, the invention proposes to use the green-to-signal ratio G _r To identify vegetation. And optimizing by using the vegetation coefficient index VI, thereby eliminating vegetation points. The green signal ratio is the ratio of the DN value of the green band to the sum of DN values of RGB; the defined formula for the vegetation coefficient VI is the result of an empirical formula. The green-to-green ratio and vegetation coefficient are defined as follows.

Wherein: VI and G _r Respectively vegetation coefficient and green letter ratio; r, G, B is the 3 color channel values of the point cloud, R ', G ', B ' are the proportions of R, G, B, and G ', B ' are defined as follows:

Notably, there are some special buildings in the real world. For example, some building roofs exist in a shed area that are blue; some residential buildings or schools have climbing tiger on the wall surface of the building, and the wall surface of the kindergarten is colored. In order to prevent false elimination of non-vegetation points during filtering of vegetation points, the invention adopts a method based on normal variation to optimize.

First, the normal vector of the point is estimated, where we use a common covariance matrix-based calculation method. If there are enough points in a region, it can construct a surface, which can be used to estimate the normal. As shown in formula (1), we set up a neighborhood N with the point p in the point cloud as the center _p . Wherein P is the original point set, P epsilon P, d (P, q) represents the distance between two points, and r is the searching radius taking P as the center of a circle.

N _p ＝{q|q∈P，d(p，q)＜r} (1)

Then we construct covariance matrix C in the neighborhood _p The definition is shown in formula (2). p is the center of all points in the neighborhood p. Calculating eigenvalues of covariance matrix and arranging lambda ₁ <λ ₂ <λ ₃ Minimum eigenvalue lambda ₁ The corresponding feature vector is considered as the normal vector to point p.

M. Pauly, point Primitives for Interactive Modeling and Processing of 3D Geometry Konstanz, germany Hartung-Gorre,2003, mentions that the maximum variance of the algorithm vector can be calculated using the eigenvalues of the covariance matrix of the normal vector, here we use this method.

We construct a neighborhood

Covariance matrix of normal vector of all points in the matrix as shown in formula (3). Same as C _p Calculating eigenvalues of the covariance matrix and arranging +.>

Can be fixedThe quantity represents the maximum variation of the normal of the point p on the gaussian sphere.

Therefore, we can set the threshold T ⁿ If at a point

A value less than T ⁿ It is considered a non-vegetation point.

1.3 building singulation

For laser point clouds, after filtering out the ground, vegetation, typically the remaining points include buildings (including fences), vehicles, and electrical facilities such as wires, poles, and the like. However, the image point cloud also includes some blocky noise points. Most non-building points can be removed by using an European clustering algorithm, and building clustering is realized.

The European clustering algorithm is a clustering algorithm based on the point spacing, and the algorithm comprises the following implementation steps:

1) Setting a threshold value T of the point spacing, wherein all points are set as non-access points;

2) The set of points P is set to be empty. Taking any one non-access point in the point cloud as a starting point, marking the point as an access point, searching all points with the space distance smaller than or equal to T from the point, and adding the points to P;

3) Repeating 2) with any one of the non-access points in the point set P as a starting point until the point set P cannot add a new point. The point set P is a cluster;

4) Repeating steps 2) -3) until all points are marked as access points.

And (5) completing clustering work.

However, some non-building points are identified as buildings, such as residual ground, planar-like green vegetation, walls connecting the buildings, etc., and the presence of these misclassifications can affect the accuracy of the roof sheet extraction. To remove these misclassifications, a Gao Chengzui small value Zmin for each class of point set is first found. In combination with the common knowledge of the roof of the building, an elevation threshold tz=zmin+4 is set. Through these steps, high-precision building singulation can be achieved.

2. Roof sheet division

After building singulation, we can take a single building as the operating unit for roof sheet extraction. Compared with an airborne laser point cloud, the building information in the image point cloud is richer, and the image point cloud not only comprises building elevation information, but also possibly comprises part of indoor information of the building.

2.1 building facade removal

Unlike on-board laser point cloud data, most buildings have complete elevation information. Common sense tells me that building facades are usually vertical to the horizontal plane. Based on this, we calculate the feature f of each point of the point cloud, defined as shown in equation (4). If the point is located in a facade, then the z value in its normal vector is most likely equal to 0, i.e. the eigenvalue f is most likely equal to 1.

f＝1-|n _p ·e _z | (4)

Because the external structure of the building in the image point cloud is clear, only the wall points of the building can be removed through the characteristic value f, and the non-wall points, such as the structures of a canopy and the like, cannot be effectively removed. The prior knowledge of the buildings tells us that the structures are distributed in small blocks and have larger distances among different blocks. Therefore, we can effectively remove by using the European clustering algorithm mentioned in 2.1.3.

2.2 division of roof sheets

Building roofs are usually made up of various planes or planes-like according to a priori knowledge, which gives us ideas. The roof sheeting may be segmented using a region growing algorithm.

The classical region growing algorithm works on the principle that the patches are segmented based on an angular comparison between point normals, i.e. by combining points that are close in smoothness by setting thresholds for curvature change and normal angle change. Thus, each output of the algorithm is considered a coplanar point. However, this approach has some drawbacks here. The algorithm will identify two parallel planes as one and the same tile, which will lead to a roof tile identification error. Meanwhile, as the method has removed the building vertical points, when the normal vector of the points is calculated, if the number of the points is considered to be more, the error of the normal vector calculation can be caused.

In combination with this information we propose a region growing algorithm incorporating the pitch of the points to segment the roof sheeting. Firstly, setting a threshold value of the point spacing, and iteratively calculating the point spacing of the nearest neighbor points, so that the point aggregation class can be used for counting the elevation mean value of each class of points. To prevent the effect of holes, we merge clusters of high Cheng Jinshi. Then we recognize different patches according to the principle of classical region growing algorithm with each class of points as input. The iteration is repeated until all points are traversed.

2.3 repair of roof sheeting

Due to calculation errors and camera shooting angles, the image point cloud often has the conditions that partial points are inconsistent with the real world, partial points are missing, and the like, for example, a transition area between a roof plane and a vertical plane is an irregular curved surface, holes are formed in a roof sheet, and the like. Meanwhile, building roofs often have some accessories such as water tanks, ventilation facilities, solar water heaters, and the like. After the vertical removal, there is still a bit of residue. These conditions will result in irregular or even incomplete boundaries of the divided roof sheets. By observing a large number of building roofs, the building roofs are usually regular panels composed of rectangles, triangles and sectors, and the points in the panels are irregularly distributed, but the boundary points have strict shape constraint. In light of this, the invention provides a new patch repair method based on PCA, which not only can regularize patches, but also can maintain the topological relation among patches.

The algorithm mainly comprises two steps: combining the roof sheet and the accessories, and repairing the sheet. The algorithm principle is as follows:

1) All point identification bits are set to 0. As shown in formula (9), the most numerous patches P are found _max The flag bit is changed to 1.

P _max ＝{P _i |max{P _i )，i∈n} (9)

Wherein: p (P) _i Represents the number of points in one panel, and n represents the number of panels on the roof.

2) Calculate all patches and P _max Is a minimum distance d of (c). To accurately find the point set of the attachment, the identification bit is selected to be 0 and 0.5<d<2.0. These patches are P _max The relevant appendages, collectively referred to as P'. Changing the identification bit of P' into 1;

3) Using principal component analysis algorithm to combine P' and P _max To the main plane. Setting the z coordinates of all points in P' as P _max Elevation mean of (c). And (5) completing projection.

4) On the principal plane, based on average dot spacing

And (5) performing point filling.

i. Find the most value X of the coordinates in the point set _min 、X _max 、X _min And X _max ；

Crossing point

Make a straight line L perpendicular to the x-axis _i Will be combined with L _i Distance is->

The points within are projected on the straight line. Calculating the distance dis between adjacent points if +.>

Then insert +.>

A point. Wherein i is 1,2,3, … up to X _min +i*d>X _max 。

ii. with (0, Y) _min ) Repeating step ii for starting point until Y _min +i*d>Y _max 。

Repeating steps ii-iii starting from the maximum of the x and y coordinates.

And (3) injection: in order not to destroy the topological relation of the patch boundary, a threshold T is set ₁ Such as T ₁ ＝4.0。

5) And restoring the P' from the main plane back to the world coordinate system to obtain the roof sheet P.

6) And iteratively executing the steps 1) to 5) until the identification bit of the point is 1.

By the processing of the algorithm, in theory, we can obtain regular and more complete building roof sheets.

2.4 refinement of roof sheeting

It has been found through experimentation that image point clouds typically contain points inside a building due to camera angles. To avoid identifying the interior points of the building as roof sheets, we propose a two-dimensional grid-based filtering method.

The principle is as follows: first, find the maximum X in the X, y axes _min 、X _max 、Y _min 、Y _max The method comprises the steps of carrying out a first treatment on the surface of the . And calculating the row and column index number of the grid to which each point belongs according to the formula (10). Counting the points in each grid, and sorting the points according to the elevation values. If the difference in elevation between adjacent points within the grid is greater than a threshold Tz, such as tz=2.5, then the portion of points having the lower elevation is discarded.

Wherein: row _i Column index representing point i _i Representing its column index, width represents the width of the grid, x _i Representing the x-coordinate, y of point i _i Representing its y-coordinate.

By adopting the technical scheme, compared with the building roof extraction technology based on the airborne LiDAR, the technical scheme of the invention has the following obvious advantages:

1) The image point cloud acquisition technology requires inexpensive, simple equipment and can obtain accurate and realistic models. The tilt imaging system mainly comprises a high-precision POS system and a plurality of cameras, and the manufacturing cost is low. In contrast, an airborne laser radar device needs hundreds of thousands of devices such as a laser scanner and an inertial navigation system, and the cost is high. Meanwhile, the image-based three-dimensional reconstruction technology can generate a precise and vivid dense point cloud model, and the extraction of the roof contour line of the building is satisfied.

2) The image point cloud data has accurate rgb information, so that the algorithm time complexity can be well reduced. Compared with LiDAR complex vegetation identification algorithm, the vegetation identification algorithm has the advantages that vegetation can be well identified through simple green-signal ratio and vegetation index calculation, so that algorithm time complexity is reduced.

3) The full-automatic building roof contour line extraction algorithm can automatically obtain the roof contour line through simple parameter setting, so that the problems of long time consumption and large investment of manpower and material resources in conventional house cadastral measurement are solved.

Drawings

FIG. 1 is a flow chart of the method of the present invention;

FIG. 2 is a schematic diagram of an original image point cloud in embodiment 1;

FIG. 3 is a schematic representation of building singulation in example 1;

FIG. 4 is a schematic view of a roof sheet point cloud in example 1;

FIG. 5 is a schematic diagram of a roof sheet point cloud in a RANSAC segmentation algorithm;

fig. 6 is a schematic diagram of a comparison of the RANSAC segmentation algorithm and the method of the present invention.

Detailed Description

Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings, wherein like or similar reference numerals refer to like or similar elements or elements having like or similar functions throughout. The embodiments described below by referring to the drawings are illustrative and intended to explain the present invention and should not be construed as limiting the invention.

In the description of the present invention, it should be understood that the terms "center", "longitudinal", "lateral", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "clockwise", "counterclockwise", etc. indicate orientations or positional relationships based on the orientations or positional relationships shown in the drawings are merely for convenience in describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, be configured and operated in a specific orientation, and thus should not be construed as limiting the present invention.

Furthermore, the terms "first," "second," and the like, are used for descriptive purposes only and are not to be construed as indicating or implying a relative importance or implicitly indicating the number of technical features indicated. Thus, a feature defining "a first" or "a second" may explicitly or implicitly include one or more such feature. In the description of the present invention, unless otherwise specified, the meaning of "a plurality" is two or more, unless otherwise clearly defined.

In the present invention, unless explicitly specified and limited otherwise, the terms "mounted," "connected," "secured," and the like are to be construed broadly and may be, for example, fixedly connected, detachably connected, or integrally connected; can be mechanically or electrically connected; can be directly connected or indirectly connected through an intermediate medium, and can be communication between two elements. The specific meaning of the above terms in the present invention can be understood by those of ordinary skill in the art according to the specific circumstances.

In the present invention, unless expressly stated or limited otherwise, a first feature "above" or "below" a second feature may include both the first and second features being in direct contact, as well as the first and second features not being in direct contact but being in contact with each other through additional features therebetween. Moreover, a first feature being "above," "over" and "on" a second feature includes the first feature being directly above and obliquely above the second feature, or simply indicating that the first feature is higher in level than the second feature. The first feature being "under", "below" and "beneath" the second feature includes the first feature being directly under and obliquely below the second feature, or simply means that the first feature is less level than the second feature.

The technical scheme of the invention is specifically described below with reference to the accompanying drawings and examples.

In order to test the correctness of the technical scheme, a group of color point cloud data generated by oblique images acquired by the unmanned aerial vehicle after a motion recovery structure algorithm is selected for experiments. The original data is collected in a certain cell of China, and total 21 txt files are different in size and comprise about 38593 ten thousand points and are about 17.6G in size. The building types in the scene are complex, and the scene comprises a kindergarten, 12 residential buildings, 2 commercial buildings with complex structures, a large number of greening facilities and vehicles, and traffic facilities such as roads and street lamps.

Because the original image has more point cloud points and larger point density, the data needs to be thinned. Meanwhile, a large number of separated block noise points exist below the ground in the scene through observation. According to the method, the thinning and clustering algorithm is executed according to file iteration, after the block noise is removed by the clustering algorithm, the thinning proportion is executed according to the ratio of 1:50, and finally the files are combined into one file, the size is 359M, and the number of points is about 772 thousands. As shown in fig. 2.

As shown in fig. 1, the implementation flow includes the following steps:

and a step a, separating the ground points from the non-ground points in the point cloud data by using a cloth filtering algorithm to obtain the non-ground points.

The specific implementation process of the algorithm is as follows:

1) Turning the original laser point cloud with outliers removed by 180 degrees;

7) Calculating the distance between the particles and the point cloud;

For specific implementation, see the relevant literature Wuming Z, jianbo Q, peng W, et al an Easy-to-Use Airborne LiDAR Data Filtering Method Based on Cloth Simulation [ J ]. Remote Sensing 2016,8 (6): 501-.

Step b, using the green-to-signal ratio G _r And identifying green or near-green points by using the vegetation coefficient index VI, optimizing the identification point set by using a variation detection algorithm based on a normal vector, and removing vegetation points.

In example step b, since vegetation is typically green or near green, the green band DN (Digital Number) values of vegetation are typically higher than the red and blue band DN values. In combination with the two points, the invention proposes to use the green-to-signal ratio G _r And vegetation coefficient index VI. The green-to-signal ratio and vegetation coefficient are defined as follows:

by setting the green-to-signal ratio G _r And a threshold value of vegetation coefficient index VI, such as G _r >0.35，VI>And 0, the green or near-green point set can be effectively identified.

However, there are often some special color buildings in the real world, such as kindergartens with near green walls in the example. In order to prevent the non-vegetation points from being removed by mistake when the vegetation points are filtered, a normal vector-based change detection algorithm is required to be adopted to further screen the identification point set, and the vegetation points are determined.

The principle of the normal vector-based change detection algorithm is as follows:

first, the normal vector of the point is estimated, here weCommon covariance matrix-based calculation methods are employed. If there are enough points in a region, it can construct a surface, which can be used to estimate the normal. As shown in formula (1), we set up a neighborhood N with the point p in the point cloud as the center _p . Wherein P is an original point set, P is P, d (P, q) represents the distance between two points, and r is a search radius taking P as a circle center.

N _p ＝{q|q∈P，d(p，q)＜r} (1)

Then we construct covariance matrix C in the neighborhood _p The definition is shown in formula (2). p is the center of all points in the neighborhood p. Calculating eigenvalues of covariance matrix and arranging lambda ₁ <λ ₂ <λ ₃ Minimum eigenvalue lambda ₁ The corresponding feature vector is considered as the normal vector of point p.

The maximum variation of the normal vector can be quantified using eigenvalues of the covariance matrix based on the normal vector mentioned in m.Pauy, point Primitives for Interactive Modeling and Processing of 3D Geometry.Konstanz,Germany:Hartung-Gorre,2003.

Specifically, construct a neighborhood

Covariance matrix of normal vector of all points in the matrix as shown in formula (3). Same as C _p Calculating eigenvalues of the covariance matrix, sorting the eigenvalues, and +. >

The maximum variation of the normal vector of the point p on the gaussian sphere can be quantitatively expressed.

By setting a threshold T ⁿ Non-vegetation points in the near-green point set can be effectively identified.

And c, carrying out connected domain analysis by using an European clustering algorithm to obtain a primary building monomalization result, and further optimizing by using a constraint term based on an elevation threshold value.

In example step c, a point cloud is obtained by step b, typically involving buildings (containing fences), vehicles and electrical practices such as wires, poles, etc. And carrying out connected domain analysis by using an European clustering algorithm, and setting a point spacing threshold and a minimum point number threshold to propose most non-building points so as to realize preliminary building clustering.

Specifically, the European clustering algorithm is a clustering algorithm based on the point spacing, and the algorithm comprises the following implementation steps:

4) Repeating steps 2) -3) until all points are marked as access points.

After the European clustering algorithm is performed, some non-building points are identified as buildings in the embodiment, such as residual ground, plane-like green vegetation, walls connected with the buildings, and the like, and the presence of the misclassifications can affect the accuracy of roof sheet extraction. To remove these misclassifications, constraint terms based on elevation thresholds need to be set here.

Specifically, a Gao Chengzui small value Zmin of each type of point set is found, and in combination with common knowledge of a building roof, an elevation threshold value, such as z > zmin+3, is set, so that most non-building point sets can be filtered. After that, the European clustering algorithm is executed again. Through these steps, high accuracy building singulation can be achieved. As shown in fig. 3 building singulation.

And executing subsequent processing steps by taking the single building point cloud as a processing unit.

And d, removing the building elevation by using the point characteristic f based on the normal vector to obtain the building roof point.

In the step d of the embodiment, in order to remove building facades of a single building, a building roof point cloud is obtained, and a point characteristic f of each point in the point cloud is calculated.

Specifically, since the building facade is generally perpendicular to the horizontal plane, the design point feature f is defined as shown in formula (4) based on this point. If the point is located on a facade, its characteristic value f is most likely equal to 1.

f＝1-|n _p ·e _z | (4)

In an embodiment, the point feature f can only remove the wall points of the building, and can not be effectively removed for non-wall points, such as structures like a canopy. And removing by adopting an European clustering algorithm.

And e, for the roof point cloud of the single building, realizing roof sheet segmentation by using a region growing algorithm based on the point spacing.

In example step e, the building roof is typically made up of various planes or planes-like, which provides a concept for us. The roof sheeting may be segmented using a region growing algorithm.

Specifically, the working principle of the classical region growing algorithm is based on angle comparison between point normals, namely, by setting a threshold value of curvature change and normal angle difference, points close in smoothness are combined, so that the patches are segmented. Thus, each output of the algorithm is considered a coplanar point. However, this approach has some drawbacks here. The algorithm will identify two parallel planes as one and the same tile, which will lead to a roof tile identification error. Meanwhile, as the method has removed the building vertical points, when the normal vector of the points is calculated, if the number of the points is considered to be more, the error of the normal vector calculation can be caused.

And f, repairing holes by using a PCA-based patch repairing algorithm on the output of the step e, and flattening the patch.

In example step f, the roof sheet obtained in step e is usually irregular or even incomplete in boundary, and even has holes. By executing the PCA-based patch repair method, the holes can be repaired, the patches can be leveled, and the topology relation between the patches can be maintained.

Specifically, the algorithm is mainly divided into two steps: combining the roof sheet and the accessories, and repairing the sheet. The algorithm principle is as follows:

P _max ＝{P _i |max{P _i )， _i ∈n) (9)

2) Calculate all patches and P _max Is a minimum distance d of (c). To accurately find the point set of the attachment, the identification bit is selected to be 0 and 0.5<d<2.0. These patches are P _max The relevant appendages, collectively referred to as P. Changing the identification bit of P into 1;

3) Using principal component analysis algorithm to sum P and P _max To the main plane. Setting the z coordinates of all points in P as P _max Elevation mean of (c). And (5) completing projection.

4) On the principal plane, based on average dot spacing

And (5) performing point filling.

i. Find the most value X of the coordinates in the point set _min 、X _max 、Y _min And Y _max ；

ii. Passing point

The points within are all projected on the straight line. Calculating the distance dis between adjacent points if +.>

Then insert +.>

A point. Wherein i is 1,2,3, … up to X _min +i*d>X _max 。

The product is prepared by (0, Y) _min ) Repeating step ii for starting point until Y _min +i*d＞Y _max 。

Repeating steps ii-iii starting from the maximum of the x and y coordinates.

Note that: in order not to destroy the topological relation of the patch boundary, a threshold T is set ₁ Such as T ₁ ＝4.0。

5) And restoring the P from the main plane back to the world coordinate system to obtain the roof sheet P.

And iteratively executing the steps 1) to 5) until the identification bit of the point is 1.

And g, carrying out de-duplication treatment on the roof sheet by utilizing the two-dimensional grid for the output of the step e.

In example step g, the output panel obtained in step f is partially a building interior panel due to camera angle and building window. For this purpose, a two-dimensional grid-based filtering method is used here to reject this partial set of points.

The principle is as follows: first, find the maximum X in the X, y axes _min 、X _max 、Y _min 、Y _max . And calculating the row and column index number of the grid to which each point belongs according to the formula (10). Statistics of each gridThe points in the table are ordered according to the elevation value. If the difference between the elevations of adjacent points in the grid is greater than a threshold Tz, such as tz=2.5, then the portion of points having the lower elevation is discarded.

The results of the roof segmentation in this example are shown in fig. 4. Experiments show that the segmentation algorithm can effectively segment the roof structure formed by rectangular combination in the data set. The analysis results were quantified by experimental data in the following manner.

Table 1 statistics of roof segmentation experimental structure

Wherein BN represents the number of the building; SP, which represents roof sheet data obtained after the roof of a building is divided; RP, which represents the number of roof panels inherent to the building in the dataset; completance, which means the integrity of the division of building roof sheets, i.e., the ratio of the number of roof sheets obtained by division to the number of roof sheets inherent to the division;

experimental results show that the roof sheet segmentation method provided by the invention can obtain better segmentation results in most buildings, has good algorithm stability and has an average segmentation result accuracy of 97.37%. Comparing the experimental results with the original data, the f and g groups of experimental results show that the under-segmentation is caused by that the two roof sheets are positioned on the same horizontal plane and are adjacent, so that the f and g groups of experimental results cannot be segmented; m and n groups of experimental data are similar in structure, and the undersegmentation occurs because the interval distance between two roof sheets is too small, so that the sheets are combined, and the number of points is too small due to the small area of the sheets, so that the segmentation is insufficient; the reason for the over-segmentation of group h is that there is a large hole in the panel and there are too many appendages nearby, thus resulting in an increased hole area after the building facade is removed, resulting in over-segmentation. These experimental results are sufficient to demonstrate the good adaptability of the algorithm of the present invention to building rooftop structures formed from convex polygon combinations.

In the specific implementation, a corresponding system is realized in a modularized mode. The method specifically comprises the following modules:

the module 1 is used for separating the ground points by using a cloth filtering algorithm to obtain non-ground points;

the module 2 is used for obtaining initial vegetation points by using the green signal ratio and the vegetation coefficient, and detecting further constraint vegetation points based on the change of the normal vector;

the module 3 is used for executing an European clustering algorithm and an elevation threshold constraint so as to realize building singulation;

a module 4 for implementing a building facade elimination algorithm and implementing roof sheet segmentation to obtain an initial roof sheet;

and a module 5, configured to implement a patch repair algorithm based on PCA and patch de-duplication based on two-dimensional grid, so as to obtain a flat roof patch.

The specific implementation of each module can participate in corresponding steps, and detailed description thereof is omitted here.

The method proposes a framework procedure involving a linear combination of algorithms. As shown in tables 2 and 3, up to 18 parameters are counted over the whole frame surface, and in practice, there are only 4 parameters to be adjusted, namely, the threshold Tf of the point feature f, the flatness threshold SmoothnessThreshold and the curvature threshold CurvatureThreshold in the region growing algorithm, the number of neighborhood points or the radius, respectively. The individual parameters are analyzed in detail below.

Table 2 parameters in building singulation techniques

Table 3 parameters in roof sheet segmentation technique

In table 2, the vegetation filter section, gr, vi is constant since the algorithm filters green plants based on the color information of the points; change threshold T for point normals ⁿ A number of experiments have found that after normalizing it based on the average value, T ⁿ =0 always achieves the desired effect, and thus can be considered as fixed. And the cloth filtering part can meet the requirements by using default parameters. Because by setting the elevation threshold T in the singulation technique _z The ground filtration results can be optimized. In the singulation, the building height is always greater than 3m, thus the elevation threshold T _z =2.5 is constant; the minimum points minP and the point spacing threshold dis are usually constant values unless the building specification or the point density varies greatly, and are minp=10000, dis=0.5, respectively.

In Table 3, the threshold T1 for the point feature f is typically about 1 for the facade removed portion, where appropriate values are required to be set by testing the sample building. For the European clustering algorithm, both the singulation and the patch segmentation are involved, the threshold of the point spacing is dis=0.5, and the minimum point number minP varies here, typically around 500. The region growing algorithm in the patch segmentation needs to properly adjust the two parameters SmoothnessThreshold, curvatureThreshold according to the quality of the image point cloud. The patch repair part Td is mainly used for preventing the damage to the topology of the building, and it generally has a minimum value, i.e., td=3, according to the common sense of the building construction; point spacing threshold T ₁ Usually of constant value, T ₁ =0.05. The patch optimizing section, the grid width is usually set to 3 times the dot pitch.

Notably, the algorithm involves multiple normal estimations, all of which employ neighborhood-based principal component analysis. Because the coordinates of the cloud point of the image point are offset compared with the real situation, the neighborhood setting must be large, and the nearest neighbor point is usually set to be 500 or the neighborhood radius is about 0.5 meter, and here, adjustment needs to be made according to the quality of the cloud point.

The method is a combined application of a series of algorithms, involving a plurality of sub-algorithms. In the prior art, only partial principle overlapping is involved, such as a region growing algorithm, an European clustering algorithm, two-dimensional meshing, a principal component analysis method and the like.

Each comparison file is analyzed in detail below to find out the prior art with true comparison meaning.

CN107220987a records a method for rapidly detecting the edge of the roof of a building based on principal component analysis, which aims to assist in detecting the edge of a roof sheet by using a principal component analysis algorithm, while the present invention uses the principal component analysis algorithm to repair the roof sheet. Is obviously different from the key points of the invention.

CN108010092B records a method for evaluating solar energy utilization potential in urban high-density areas based on low-altitude photogrammetry data modeling. In this document, a building singulation section attempts to separate a building point cloud from color information in the image point cloud by means of volvox plug-ins in a parameterized modeling tool grasshoper. It is disadvantageous in that it does not take into account that the color of the building is similar to that of the ground or vegetation. Is obviously different from the key points of the invention.

CN107944384a (covariance matrix-region growing) discloses a building roof tile segmentation method based on a three-dimensional Voronoi diagram, and the processing object is an airborne LiDAR point cloud. According to the method, the three-dimensional Voronoi diagram is used for replacing kdtree to establish the point neighborhood, so that parameter setting during neighborhood setting is avoided, and the problem that the neighborhood is difficult to control during construction of the point cloud data space topological relation is effectively solved. However, as the processing object is an airborne LiDAR point cloud, the vegetation filtering algorithm adopted by the processing object is completely different from the invention, and the processing object is obviously different from the core key point of the invention.

CN107230251a discloses a technique to create a 3D city model from oblique imaging data and lidar data. This document relates to a method of matching and fusing oblique imaging data and lidar data to create a 3D city model and design of a hybrid 3D imaging device. The invention describes an algorithm for processing a color three-dimensional point cloud, and the algorithm are not contrasted. The key points of the invention are obviously different from those of the invention.

CN108090957a discloses a method of mapping terrain based on BIM. The document describes a method for establishing a three-dimensional point cloud model mainly based on man-machine interaction, and the invention relates to a full-automatic building roof sheet recognition algorithm which is obviously different from the key points of the invention.

CN106846494a (grid method_gpu) discloses an automatic monomer algorithm for oblique photography three-dimensional building models. The algorithm enables building singulation by passing vertex information of the building contours into the GPU shader, and by highlighting the building markers within the contours. This document uses GPU technology, whose technical route is very different from that of the present invention. Is obviously different from the key points of the invention.

CN108074232a (grid-element) is the same as CN108109139a, and an airborne LiDAR building detection method based on element segmentation is disclosed. The method comprises the steps of regularizing original airborne LiDAR point cloud data into a gray 3D volume metadata set, and dividing and marking the gray 3D volume metadata into a plurality of 3D connected areas based on connectivity and radiation characteristic similarity criteria; then, based on the characteristics of the roof and the elevation of the building, the detection of the airborne LiDAR building based on voxel segmentation is completed. The method adopts a voxel technology and is obviously different from the key point of the invention.

CN104809689B discloses a contour-based map registration method for building point cloud models. The aim of this document is to achieve registration between the point cloud model and the satellite base map. Is obviously different from the key points of the invention.

CN102411778B discloses an automatic matching method of airborne laser point cloud and aerial image. The method directly extracts a building contour line from an airborne LiDAR point cloud, acquires building angle features of registration primitives based on the building contour line, and then automatically matches homonymous angle features between the point cloud and an aerial image with the assistance of an external azimuth element; and then, adopting a beam method area network adjustment and loop iteration strategy to realize the overall optimal registration of the aerial image and the point cloud data. Is obviously different from the key points of the invention.

CN101726255B discloses a method of extracting a building of interest from three-dimensional laser point cloud data. This approach attempts to separate buildings from non-ground points directly using the European clustering algorithm, which has certain limitations in theory. Because non-building points such as vegetation, large automobiles and the like are often arranged on the ground in a real scene, the robustness of simply executing the European clustering algorithm is not high. At the same time, the object of this document is to search for buildings by setting building edge features. Thus, it is clearly different from the core point of the present invention.

CN107545602a (spatial topological relation-building modeling) discloses a building modeling method under spatial topological relation constraint based on LiDAR point cloud. This document focuses on the spatial topological relation handling of building roof geometry primitives with little detail in building singulation and roof segmentation. Is obviously different from the key points of the invention.

CN104036544B discloses a building roof reconstruction method based on airborne LiDAR data. The key point of the document is how to accurately obtain the vector boundary of the roof sheet, the segmentation part of the roof sheet of the building is simpler, the roof sheet is obtained by directly adopting a region growing algorithm, and then a plane equation of the plane where the roof sheet is located is obtained by utilizing least square fitting. Is obviously different from the key points of the invention.

CN105572687B discloses a method for making a digital line drawing of a building based on a vehicle-mounted laser radar point cloud. The method is used for detecting and extracting the elevation of the building by the vehicle-mounted laser radar data, and has no contrast with the method.

The tree extraction method based on region growing and gradient segmentation is not compared with the method.

CN106970375a discloses a method for automatically extracting building information in an airborne laser radar. The method essentially utilizes a plane fitting method to detect all the building point cloud planes, and then utilizes three-dimensional morphological corrosion operation to remove part of misclassification point clouds, so that accurate building point clouds are obtained. The method is applicable to airborne LiDAR point clouds, however, unlike laser point clouds, the image point cloud is high in point cloud density, and contains a complete building elevation point while containing a roof point cloud. Therefore, if no other measures are taken, the plane fitting method is directly performed, and the building point cloud cannot be detected through a small number of iterations. Thus, it is clearly different from the core point of the present invention.

A building roof cloud plane segmentation method based on local constraint introduces a RANSAC algorithm added with point cloud point normal vector constraint, and solves the problem of separation of different roof surfaces on the same plane of the traditional RANSAC algorithm to a certain extent. The patch segmentation algorithm adopts a region growing algorithm, and is obviously different from the key points of the invention.

The journal paper (improved RANSAC point cloud segmentation algorithm considering the roof structure of a building) utilizes a triangle region growth method of gradient and height difference to decompose different structural layers of a complex building, and then provides a RANSAC algorithm with a floating consistent set threshold value to detect the roof plane of the building, so that the method has certain applicability. The patch segmentation algorithm adopts a region growing algorithm, and is obviously different from the key points of the invention.

The research graduate paper on-board LiDAR point cloud data filtering and building point group segmentation research improves the extraction method of random Hough transformation on the point cloud contained in the building roof surface patch. The assumption of the whole algorithm is based on that the building model is composed of a plurality of planes, which has a certain limitation, because in a real scene, some of the buildings are often composed of curved surfaces. Is obviously different from the key points of the invention.

CN109242855a discloses roof segmentation method, system and device based on multi-resolution three-dimensional statistical information. According to the method, three-dimensional point cloud feature statistical information with different resolutions is extracted from the image point cloud, and then semantic classification is carried out on a three-dimensional point cloud scene by directly utilizing global energy optimization, so that the point cloud of a building is obtained. The whole technical route is different from the method adopted by the invention.

CN108171720a discloses a method for detecting the boundary of an object of an oblique photography model based on geometric statistics. According to the method, a two-dimensional minimum outsourcing rectangle of the building is calculated, the minimum outsourcing rectangle is intersected with a model, and point cloud data belonging to the area where the current object is located is selected to obtain the individualized building. The method has complex process and large calculation amount. Is quite different from the technical route adopted by the invention.

CN108898144a discloses a building damage status detection method. The method directly takes building point cloud data as an input object, and does not relate to building segmentation, roof detection and the like, and has no contrast.

CN109461207a discloses a method and device for the singulation of point cloud data buildings. Before the voxels are classified, the method carries out pretreatment on three-dimensional points in the voxels, and the three-dimensional points are divided again, so that the three-dimensional points contained in each voxel are more similar. The method uses a global energy optimization function, and the calculated amount is large. Is quite different from the technical route adopted by the invention.

CN109754020a discloses a ground point cloud extraction method integrating a multi-level progressive strategy and unsupervised learning. The method is a ground point cloud extraction algorithm, and is obviously different from the key points of the invention.

CN109859315a discloses a method for separating vegetation on the ground surface of a three-dimensional image. The method is a software operation method and is not contrasted with the present invention.

CN109870106a discloses a building volume measuring method based on unmanned aerial vehicle pictures. The method comprises the steps of obtaining building image point clouds through intensive matching of unmanned aerial vehicle images, and then calculating the volume of a building by an integration method for a model after Delaunay triangulation. Is obviously different from the key points of the invention.

CN105844629B (RGA-DP-RANSAC) discloses an automatic segmentation method for the standing point cloud of large-scene city buildings. The method for extracting the cloud data of the roof point of the airborne LiDAR building comprises the following steps: and (3) realizing the separation of ground points and non-ground points by adopting progressive irregular triangular network encryption, filtering ground feature points with the height difference smaller than 2.0m by taking the ground point height as a reference, and then dividing building roof point cloud by adopting a RANSAC (random sample area network) face patch detection algorithm. The method is relatively simple and is not applicable to image point clouds with a large amount of elevation information. Is obviously different from the key points of the invention.

CN106097311a (progressive morphological filtering-region growing-minimum loop detection) discloses a building three-dimensional reconstruction method of airborne laser radar data. In the building point cloud extraction, the number of point cloud echoes is not applicable here because it is not found in the image point cloud. Second, it employs a conventional region growing algorithm to identify building point clouds, which has certain limitations. And (3) in the building roof segmentation part, after a roof plane is detected by adopting a clustering growth segmentation algorithm and a RANSAC algorithm based on the point cloud space distribution characteristic, the optimization is carried out by utilizing the panel normal vector included angle constraint and the panel distance constraint, and similar panels are combined. The key points of the invention are obviously different from those of the invention.

CN106600680a discloses a batch fine three-dimensional modeling method for building object frame models. The method is a set of mature production and quality inspection operation flow. Is obviously different from the key points of the invention.

CN107644452a (airborne LiDAR-roof tile segmentation) discloses an airborne LiDAR point cloud roof tile segmentation method and system. The method is different from the invention in that:

1) And constructing a point cloud neighborhood system by taking the 3D Voronoi diagram as airborne LiDAR point cloud data, wherein the method comprises the step of constructing an adjacency relationship based on the 3D Voronoi neighborhood for any point in the airborne LiDAR point cloud data. The present invention employs kdtree.

2) This document utilizes multiple returns of the LiDAR point cloud to distinguish building points from vegetation points. The processing data of the invention is image point cloud, no echo information is generated, and other methods are adopted.

CN104484668B discloses a method for extracting building contour lines of multi-overlapping remote sensing images of unmanned aerial vehicle. The method is different from the invention in that:

1) The vegetation in the ground points is filtered by adopting the color invariant, and the situation that the colors of the buildings are similar is not considered. After the color filtering is adopted, the method adds a normal vector-based change detection algorithm, so that the color filtering is ensured not to cause the loss of building points.

2) A building point cloud is obtained using a planar fit-based region growing method, which is essentially based on the assumption that the building is made up of planes. After vegetation and the ground are filtered, the European clustering algorithm is directly adopted, and the building is obtained by setting a point spacing threshold value and the minimum clustering point number. It is not based on this assumption, ensuring that buildings of arbitrary structure can be identified.

3) The facade splitting section uses the direction of the normal vector of the face sheet to split the building facade. The invention uses normal vectors of points.

CN106023312B (facade removal-plane fitting-region growing) discloses an automatic reconstruction method for three-dimensional building model based on aviation LiDAR data. The method is different from the invention in that:

1) The method that the density of building elevation points in the airborne LiDAR point cloud is smaller than that of roof points is grasped to provide wall points;

2) The roof layer resampling part uses a support vector machine algorithm, so that the calculated amount is increased.

3) The seed area selection-roof sheet growth-sheet leveling optimization part adopts the method that initial plane parameters are preliminarily obtained by utilizing a point curvature threshold value, and then the roof sheet growth is realized by utilizing a distance threshold value and a distance standard deviation.

In summary, the present invention differs from the prior art in that:

in building identification, the elimination method is employed. Namely, the ground object is divided into 3 types: building, ground and vegetation, removing ground points and vegetation points, and the rest is building points. Many of the methods in the references are based on the assumption that a building is made up of planar panels, which has certain limitations that do not allow for extraction of buildings that are not planar. The present invention does not suffer from this problem.

In the aspect of roof sheet segmentation, after the original roof sheet point cloud is extracted, a face sheet repairing method based on PCA is adopted, so that the problem of holes in the original image point cloud is solved, and the face sheet is flattened. In the reference document, the processing of holes in the image point cloud is not involved.

In particular, the method comprises the steps of,

1) In the ground filtering algorithm, the latest cloth filtering algorithm is adopted. The method has the advantages of stable algorithm, simple parameter setting, easy understanding of principle and the like, and meanwhile, the method has strong scene applicability, can treat flat ground and is also applicable to the ground with a certain gradient. While most of the references use progressive morphological filtering algorithms.

2) In the vegetation filtering algorithm, a new method is invented. First, the point cloud is divided into green and non-green points based on the color information of the points. Then, the variation of the normal vector in the green point measured by the eigenvalue in the covariance matrix is calculated to determine which points belong to the building point. Under the dual actions of the front step and the rear step, green vegetation can be removed, and meanwhile, the building points can not be deleted by mistake. There are also methods in the literature that use color information to filter vegetation similarly, but they do not consider the situation where building colors are similar to vegetation, and thus algorithms have certain limitations.

3) In a roof patch repair algorithm, a patch repair method based on PCA projection conversion is provided, which can process accessories of a roof, make patches smooth, repair holes and maintain topological relations among patches. The problem of repairing holes in the patch is not considered in the reference.

4) In the refinement of the roof sheet, a method for identifying overlapping sheets based on a two-dimensional grid is adopted. In airborne LiDAR point cloud or unmanned aerial vehicle inclined image matching point cloud, because of the inclination angle problem of a laser radar or a camera, the point cloud often contains some building interior points and even a building interior panel point set, and here, the deletion of the building interior points is completed by eliminating the mode of the lower elevation of the grid point set by carrying out the high Cheng Ju class on the grid point set, so that the accuracy of identifying the building roof panel is improved.

The comparison of the inventive method with the prior art examples is verified below.

1. The aspect of vegetation filtering algorithm:

the color invariant based plant filtering method mentioned in CN104484668B was chosen here for comparison with the method of the present invention and the test dataset was the dataset in the specific example.

Specifically, let the coordinates of each point in the point cloud be (x, y, z), the three color channels be (R, G, B), the threshold for vegetation of the color invariant be Tg, and the color invariant formula defined by the green and blue channels be:

wherein Ig (x, y, z), I _b (x, y, z) represents the green and blue component values of the point cloud at the (x, y, z) point. Psi phi type _g (x, y, z) represents the color invariant at the (x, y, z) point. When psi is _g ＜T _g When the point is indicated as a vegetation point; otherwise, it is a non-vegetation point.

Table 3 vegetation filtering algorithm comparison

Method	Pts
		The method of the invention	216335
Color invariant theory	225893

Where Pts represents the vegetation point obtained by the filtering algorithm. Since the number of correct vegetation points cannot be quantified, we can only compare results according to visual results of naked eyes. The comparison experiment shows that under the condition that the correct vegetation points identified by naked eye observation are the same, the method adds a normal vector-based change detection algorithm, so that the near-green building point set is well protected, and the number of the result points of the method is smaller than the color invariance theory. In conclusion, the algorithm of the invention has better stability and is due to the existing vegetation filtering algorithm based on the color information.

2. Roof sheet segmentation algorithm aspect:

here we choose the RANSAC point cloud segmentation algorithm to compare with the algorithm proposed in this embodiment. The segmentation effect is shown in fig. 5. The results were as follows:

table 4 statistics of roof segmentation experimental structure

For better comparison of the two algorithms, three buildings g, m and n were chosen for individual analysis. As shown in fig. 6, the roof sheet extraction results of three buildings g, m, and n are shown. Wherein (a-c) is the RANSAC method and (d-f) is the method of the invention. The result shows that the RANSAC point cloud segmentation algorithm is sensitive to the number of the surface patches, and the small surface patches cannot be effectively extracted. Therefore, the algorithm of the invention has better stability and can better keep the detail structure of the roof surface.

In the description of the present specification, a description of the terms "one embodiment," "some embodiments," "examples," "specific examples," or "some examples," etc., means that a particular feature, structure, material, or characteristic described in connection with the embodiment or example is included in at least one embodiment or example of the invention. In this specification, schematic representations of the above terms do not necessarily refer to the same embodiments or examples. Furthermore, the particular features, structures, materials, or characteristics described may be combined in any suitable manner in any one or more embodiments or examples.

Although embodiments of the present invention have been shown and described above, it will be understood that the above embodiments are illustrative and not to be construed as limiting the invention, and that variations, modifications, alternatives, and variations may be made in the above embodiments by those skilled in the art without departing from the spirit and principles of the invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for partitioning a building roof sheeting from a large-scale image dense matching point cloud, comprising the steps of:

3) Removing the building elevation;

4) Cutting, repairing and refining the roof sheet;

the division of the roof sheet refers to the division of the roof sheet by utilizing a region growing algorithm based on the point spacing for the roof point cloud of a single building; the repairing of the roof sheet refers to the steps of repairing holes by using a PCA-based sheet repairing algorithm on output of the segmented roof sheet and flattening the sheet; the refinement of the roof sheet refers to the output of the divided roof sheet, and the two-dimensional grid is utilized to realize the de-duplication treatment of the roof sheet;

The PCA-based patch repair algorithm comprises the following specific steps:

1) Setting all the point identification bits to 0; as shown in the following formula (9), the most numerous patches P are found _max The identification bit is changed into 1;

P _max ＝{P _i |max{P _i }，i∈n} (9)

wherein: p (P) _i Represents the number of points in one panel, n represents the number of panels on the roof;

2) Calculate all patches and P _max A minimum distance d of (2); to accurately find the point set of the attachment, the identification bit is selected to be 0 and 0.5<d<2.0 facesA sheet; these patches are P _max Related appendages, collectively referred to as P'; changing the identification bit of P' into 1;

3) Using principal component analysis algorithm to combine P' and P _max Switching to a main plane; setting the z coordinates of all points in P' as P _max Elevation means of (2); completing projection;

4) On the principal plane, based on average dot spacing

Performing point compensation;

5) Restoring P' from the main plane back to the world coordinate system to obtain a roof sheet P;

6) And iteratively executing the steps 1) to 5) of the patch repair algorithm based on the PCA until the identification bit of the point is 1.

2. The method of claim 1, wherein the ground filtering separates ground points from non-ground points in the point cloud data using a cloth filtering algorithm to obtain non-ground points.

3. The method for partitioning building roof tiles from a large-scale image dense matching point cloud of claim 1, wherein the ground filtering is performed as follows:

1) Turning the original laser point cloud with outliers removed by 180 degrees;

3) Projecting the laser point cloud and the grid particles to a horizontal plane, finding the nearest neighbor laser point of each particle, and taking Gao Chengji of the laser points as an intersection point elevation value;

7) Calculating the distance between the particles and the point cloud;

4. A method of segmenting building roof tiles from a large scale image dense matching point cloud as claimed in claim 3, wherein the ground filtering is performed in a process further comprising a post-processing step of: searching for the four neighboring particles of each movable particle, and if an immovable particle is found, comparing the two to the elevation of the corresponding spot; if the minimum elevation difference is less than the threshold value, the movable particles continue to move until set to be immovable; this step is iteratively performed until all movable particles are traversed.

5. The method of segmenting building roof tiles from a large scale image dense matching point cloud as recited in claim 1, wherein the vegetation filtering utilizes a green to blue ratio G _r And vegetation coefficient index VI to identify green or near-green points, and optimizing the identification point set by utilizing a variation detection algorithm based on normal vectors to further remove vegetation points _。

6. The method of claim 5, wherein the green signal ratio is a ratio of DN values of green bands to a sum of DN values of RGB; the green-to-signal ratio and vegetation coefficient are defined as follows:

7. the method for partitioning building roof tiles from a large-scale image dense matching point cloud as recited in claim 1, wherein said building is singulated, a connected domain analysis is performed by using an euclidean clustering algorithm to obtain a preliminary building singulation result, and then a constraint term based on an elevation threshold is used for further optimization.

8. The method for partitioning building roof tiles from a large-scale image dense matching point cloud of claim 1, wherein removing building facades uses normal vector based point feature f to remove building facades to obtain building roof point clouds; the definition of the point feature f is shown in the following formula:

f＝1-|n _p ·e _z |

9. A system for segmenting building roof sheeting from a large scale image dense matching point cloud, employing the method of any of claims 1 to 8, comprising: