CN110570428A

CN110570428A - method and system for segmenting roof surface patch of building from large-scale image dense matching point cloud

Info

Publication number: CN110570428A
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: 2019-12-13
Anticipated expiration: 2039-08-09
Also published as: CN110570428B

Abstract

The invention discloses a method and a system for dividing a building roof surface patch from a large-scale image dense matching point cloud, which are used for emphatically solving the problem of accurately acquiring the building roof three-dimensional surface patch from a color three-dimensional point cloud generated from an oblique image. A method for segmenting a roof surface patch of a building from a large-scale image dense matching point cloud comprises the following steps: ground filtration and vegetation filtration; removing the ground and vegetation in the image point cloud; the building is integrated; eliminating non-building points and realizing building clustering; removing the building facade; and (4) dividing, repairing and refining the roof surface patch.

Description

Method and system for segmenting roof surface patch of building 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 segmenting a roof surface patch of a building from a large-scale image dense matching point cloud.

Background

cities are changing day by day as important subjects of spatio-temporal information. The three-dimensional model of the urban building is used as a visual expression of information in a series of urban spaces such as urban landforms, buildings on the ground and the like, shows various information such as geometric forms, attributes, spatial positions, textures and the like of a target, and is the basis of an urban three-dimensional geographic information system. The method has important functions in various aspects such as GIS database updating, land utilization, digital cities and the like by quickly extracting the buildings and detecting the change of the buildings.

The method takes three-dimensional color point cloud obtained by an oblique photography measurement system as a data source to research the roof surface segmentation technology of the building. Details of the oblique photogrammetry system and the building rooftop extraction technology based on airborne LiDAR are described below.

The oblique photogrammetry system is a new technology developed on the basis of the traditional photogrammetry technology by combining the oblique photogrammetry technology of computer vision, and can simultaneously acquire ground images of a plurality of angles and acquire images of the top and rich textures of the side surfaces of a building. The technology has the advantages that the real situation and form of the ground object can be presented in an all-around and all-angle manner, the high-precision side texture information of the ground object, particularly a building, can be captured, the ground object in the image can be positioned and modeled, a real urban three-dimensional model is constructed, and the application field of the technology becomes more and more extensive.

at present, the tilt imaging system is basically composed of a "1 + N" camera, i.e. one downward-looking camera and N tilt cameras, and the "1 + N" tilt imaging system can acquire images with complete coverage from the perspective of urban three-dimensional modeling analysis. Currently, 4 oblique cameras, i.e., 2 cameras in front and back, left and right, are most common. The inclined imaging system is basically integrated with a high-precision POS system, and is used for acquiring images and recording the position and posture data of a camera at the time of image exposure.

in order to further improve the completeness of the image coverage of the tilt imaging system, the adopted technical strategy is to increase the number of tilt cameras, such as the octobluque MIDAS of the Track Air formula, to 8, and to shoot in a 360-degree panorama by one downward-looking camera and eight tilt cameras with a tilt angle of 45 degrees, and the additional four cameras can create favorable conditions for covering the shot area in an all-round and non-dead-angle manner.

In summary, the multi-view oblique camera is generally mounted on an aviation aircraft to acquire ground feature images from multiple angles. From the technical characteristics of multi-view oblique photography, the method mainly has the following 4 characteristics:

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

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

(3) The vertical image and the rear-view image of the same route with an overlapping area have a rotation of about 180 degrees;

(4) In the absence of a cross-course, there is a rotation of approximately ± 90 ° between the vertical image and the left-right view image with the overlapping area.

the technology for three-dimensional reconstruction of dense point cloud based on image, also called multi-view stereo (MVS), uses a calibration image obtained from SFM algorithm, i.e. parameters of a camera corresponding to an image are obtained by calculation or other methods. At present, the post-processing business software mainly comprises ContextCapture software produced by Acute3D, France, an Altizure platform established by the professor of hong Kong science and technology university, and the like. The three-dimensional reconstruction technology of the dense point cloud based on the image comprises the following main processes:

(1) detecting and matching the feature points;

(2) Solving camera parameters; usually, a camera solving method in the SFM meaning is used, and under the inaccurate initialization condition, bundle adjustment constraint is utilized to obtain accurate camera parameters and obtain a sparse point cloud model;

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

(4) Mesh hole patching and adding texture mapping.

The oblique photogrammetry technology combined with computer vision breaks through the limitations of the remote sensing technology, the traditional photogrammetry and surveying and mapping industries, 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 'smart cities' and 'digital cities'. In the method, the oblique photography technology is the basis of research, and the dense color three-dimensional point cloud obtained by the method is the object of the extraction of the roof contour line of the building.

Generally speaking, building rooftop extraction techniques based on airborne LiDAR are primarily 2 steps: building singlelization and roofing 3D patch generation. Two approaches are commonly used for building singleness. 1) Identifying buildings directly from the original point cloud; 2) and respectively filtering the ground and the vegetation to obtain the building point cloud. The roof surface patch extraction technology is carried out according to building facade removal and surface patch segmentation.

disclosure of Invention

In order to solve the above technical problems, an object of the present invention is to provide a method and a system for segmenting a building roof surface patch from a large-scale image dense matching point cloud, and to emphasize how to accurately obtain a building roof three-dimensional surface patch from a color three-dimensional point cloud generated from an oblique image.

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

A method for segmenting a roof surface patch of a building from a large-scale image dense matching point cloud comprises the following steps:

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

2) the building is integrated; eliminating non-building points and realizing building clustering;

3) Removing the building facade;

4) And (4) dividing, repairing and refining the roof surface patch.

Preferably, the ground filtering is to separate ground points from non-ground points in the point cloud data by using a cloth filtering algorithm to obtain the non-ground points.

preferably, the ground filtration is realized by the following specific steps:

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

2) Initializing a cloth grid; determining the total number of the 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 recording the elevation of the laser point as an intersection point elevation value;

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

5) For each mesh particle, calculating a positional shift due to interaction forces between the particles;

6) Iteratively performing the above steps 4) to 5) 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 particle and the point cloud;

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

Preferably, the implementation process of the ground filtration further comprises a post-processing step: searching four adjacent particles of each movable particle, and if finding that the immovable particles exist, comparing the elevations of the movable particles and the immovable particles to the stress light spot; if the minimum elevation difference is less than the threshold, the movable particle continues to move until set to be immovable; this step is performed iteratively until all movable particles are traversed.

Preferably, the vegetation filtration utilizes the split G_rIdentifying green or near-green color points according to the vegetation coefficient index VI, and then utilizing a change detection algorithm based on a normal vectorAnd optimizing the identification point set, and further removing the vegetation points.

Preferably, the green signal ratio is the ratio of DN value of green wave band to the sum of DN values of RGB; the definition of the split and vegetation coefficient is shown as follows:

in the formula: VI and G_rRespectively vegetation coefficient and split green; r, G, B are the 3 color channel values of the point cloud, R ', G', B 'are the proportion of R, G, B, R', G ', B' are defined as the following formula:

Preferably, the building is individualized, a connected domain analysis is performed by using an Euclidean clustering algorithm to obtain a preliminary building individualization result, and then a constraint term based on an elevation threshold is used for further optimization.

Preferably, removing the facade of the building, namely removing the facade of the building by using the point feature f based on the normal vector to obtain the point cloud of the roof of the building; the definition of the point feature f is shown as follows:

f＝1-|n_p·e_z|

In the formula, n_pIs the normal vector at point p, e_zis the unit column vector (0,0, 1).

Preferably, the roof surface patch is divided by utilizing a region growing algorithm based on point spacing to the single building roof point cloud; the repairing of the roof surface patch refers to that the hole is repaired by utilizing a patch repairing algorithm based on PCA (principal component analysis) on the output of the divided roof surface patch, and the patch is flattened; the refinement of the roof surface patch refers to that the two-dimensional grid is utilized to realize the duplicate removal processing of the roof surface patch for the output of the divided roof surface patch.

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

The first module is used for separating 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 coefficients and detecting further constrained vegetation points based on the change of the normal vector;

The third module is used for executing an Euclidean clustering algorithm and elevation threshold value constraint so as to realize building singleness;

The fourth module is used for realizing a building facade elimination algorithm and realizing roof surface patch segmentation to obtain an initial roof surface patch;

And the fifth module is used for realizing a patch repairing algorithm based on PCA and patch duplication removal based on the two-dimensional grid to obtain a flat roof patch.

By using some ideas of LiDAR point cloud processing, the invention provides an algorithm framework for fully automatically extracting three-dimensional contour lines of buildings from image point clouds. Because the image point cloud has a lot of noises, compared with the LiDAR point cloud, the algorithm is more complex. According to the flow, the framework mainly comprises three parts: building singleization, roof surface slice 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 singleton

In order to accurately extract the three-dimensional contour line of the building from the large-scale image point cloud, the building is required to be integrated. Because the roof contour line is extracted by taking a single building as a unit, the contour line is extracted on the basis of the unit. By means of the idea of 'bottom-up' separation of buildings by LiDAR point cloud, the method realizes the building monomer formation by sequentially separating ground and vegetation.

1.1 ground filtration

the image point cloud is similar to the airborne laser point cloud by analyzing the data properties, so that the ground can be removed by using a ground identification algorithm in the airborne laser point cloud. Here, the Cloth filter algorithm proposed by Wuming Zhang in 2016 is used, and the paper is detailed in Wuming Z, Jianbo Q, Peng W, et al, an Easy-to-Use air filtration LiDAR data filtration Method Based on Cloth filtration [ J ]. removal Sensing,2016,8(6): 501-. The algorithm has the advantages of less parameter setting, strong algorithm generalization capability and capability of obtaining stable ground filtering effect.

The cloth filter algorithm simulates a simple physical process. Suppose a piece of cloth falls in the air due to gravity. If the cloth is soft enough to adhere to the surface of an object, the final shape of the cloth is a digital surface model of the scene. The point cloud scene is now flipped 180 degrees so that the ground portion in the point cloud is above. In contrast, here the cloth has a certain stiffness, and the resulting shape of the cloth is a digital terrain model, i.e. the ground. Analysis of the interaction of the particles of the cloth and the corresponding laser point cloud allows the final shape of the cloth to be determined, which is used to separate ground points from non-ground points. Based on the idea, the algorithm uses a cloth simulation technology in three-dimensional computer graphics, namely, cloth is simulated through computer programming. The cloth here is a grid formed by a mass of particles of a fixed mass but no spatial size connected. The vertices of the mesh are these physical particles and the edges of the mesh are a "spring" that obeys Hooke's Law.

the specific implementation process of the algorithm is as follows:

2) And (5) initializing a cloth grid. Determining the total number of the cloth particles by the grid resolution defined by a user; the initial position of the cloth is placed above the highest point of the point cloud;

3) and projecting the laser point cloud and the grid particles to a horizontal plane. Finding the nearest neighbor laser point of each particle, and recording the elevation of the laser point as the elevation value of the intersection point;

4) for each particle, its position due to gravity is calculated. If its elevation is greater than the intersection elevation, it may continue to move; if the height value of the particle is less than or equal to the intersection point height value, the height value of the particle is set as the intersection point height value and marked as immovable;

6) iteration is performed 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 particle and the point cloud;

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

in order to adapt the algorithm to terrain with large gradient changes, a post-processing algorithm is added. The principle is as follows: searching four adjacent particles of each movable particle, and if the immovable particles are found, comparing the elevations of the two on the stress light spot. If the minimum elevation difference is less than the threshold, the movable particles continue to move until set to be immovable. This step is performed iteratively until all movable particles are traversed.

1.2 Vegetation filtration

The algorithm does not adopt the idea of traditional laser point cloud. From data content analysis, the data source is a color three-dimensional point cloud with correct RGB information; since vegetation is typically green or near-green, the DN (DigitalNumber) values of the green band of vegetation are generally higher than the DN values of the red and blue bands. In combination with these two points, the present invention proposes the use of the split G_rTo identify vegetation. And optimizing by utilizing the vegetation coefficient index VI, thereby realizing the elimination of vegetation points. The green signal ratio is the ratio of DN value of green wave band to the sum of DN values of RGB; the formula for defining the vegetation coefficient VI is the result of an empirical formula. The definition of the split and vegetation coefficient is shown below.

In the formula: VI and G_rRespectively vegetation coefficient and split green; r, G, B are the 3 color channel values of the point cloud, R ', G ', B ' are the proportion of R, G, B, G ', B ' are defined as the following formula:

It is worth noting that there are some special buildings in the real world. For example, there are some buildings in the shed area that have roofs that are blue; some residential buildings or schools have mountain climbing tigers on the wall surfaces, and the wall surfaces of kindergarten are colored. In order to prevent non-vegetation points from being mistakenly removed when the vegetation points are filtered, the method adopts a method based on normal line change to carry out optimization.

First, the normal vectors of the points are estimated, where we use a common covariance matrix based calculation method. If there are enough points in a region that it can construct a surface, we can use it to estimate the normal. As shown in formula (1), a neighborhood N is established by taking a point p in a point cloud as a center_p. Wherein, P is an original point set, P belongs to P, d (P, q) represents the distance between two points, and r is a search radius taking P as the center of a circle.

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

Then, we construct the covariance matrix C in the neighborhood_pThe definition is shown in formula (2). p is the center of all points in the neighborhood p. Calculating eigenvalues of the covariance matrix, and ordering λ₁<λ₂<λ₃Minimum eigenvalue λ₁the corresponding feature vector is considered to be the normal vector for point p.

Pauly, Point principles for Interactive Modeling and Processing of 3D geometry, Konstanz, Germany: Hartung-Gorre,2003, mentions that the maximum variation of the algorithm vector can be calculated using the eigenvalues of the covariance matrix of the normal vector, where we use this method.

We construct neighborhoodscovariance matrix of normal vectors of all points insidematrix, as shown in equation (3). Same as C_pcalculating eigenvalues of the covariance matrix and orderingthe maximum variation of the normal of the point p on the gaussian sphere can be quantified.

therefore, we can set the threshold TⁿIf appropriate, ofvalue less than TⁿIt is considered as a non-vegetation point.

1.3 building singlets

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

the Euclidean clustering algorithm is a clustering algorithm based on point spacing, and the algorithm implementation steps are as follows:

1) Setting a threshold T of dot spacing, all dots being set as unaccessed dots;

2) set point set P to null. 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 a space distance less than or equal to T from the point, and adding the points to P;

3) Repeat 2) starting from any one of the unaccessed points in the point set P until the point set P cannot be added with new points. The point set P is a cluster;

4) Repeat execution 2) -3) until all points are marked as access points.

And finishing clustering work.

however, some non-building points are identified as buildings, such as residual ground, plane-like green vegetation, walls connecting buildings, and the like, and the presence of such misclassification affects the accuracy of roof tile extraction. To remove these misclassifications, the minimum elevation value Zmin of each point set is first found. And setting the elevation threshold value Tz to Zmin +4 by combining the common knowledge of the building roof. Through these steps, the building can be integrated with high precision.

2. Roof panel segmentation

after the building is formed into a single body, the roof surface slice can be extracted by taking a single building as an operation unit. Compared with airborne laser point cloud, the building information in the image point cloud is richer, and the image point cloud not only comprises building facade information, but also possibly comprises part of indoor information of the building.

2.1 building facade removal

The image point cloud data is different from the airborne laser point cloud data, and most buildings have complete facade information. Conventional wisdom tells us that the building facade is usually vertical to the horizontal. Based on this, we calculate the feature f of each point in the point cloud, defined as formula (4). If a point is located on a facade, it is very likely that the value of z in its normal vector will equal 0, i.e., the eigenvalue f will be very likely equal to 1.

f＝1-|n_p·e_z| (4)

because the exterior surface structure of the building in the image point cloud is relatively clear, the wall surface points of the building can only be removed through the characteristic value f, and the non-wall surface points, such as structures like a canopy and the like, can not be effectively removed. The a priori knowledge of buildings tells us that these structures are distributed sporadically in small blocks, with large distances between different blocks. Therefore, we can effectively remove it by using the Euclidean clustering algorithm mentioned in 2.1.3.

2.2 splitting of roof tiles

According to the prior knowledge, the building roof is generally composed of various planes or plane-like planes, which provides a thought for us. The roof patch may be segmented using a region growing algorithm.

The working principle of the classical region growing algorithm is based on the angular comparison between point normals, i.e. by setting thresholds for curvature change and normal angle change, merging points that are close in smoothness, thus dividing a patch. Thus, each time the algorithm outputs, it is considered to be a co-planar point. However, this approach has some disadvantages here. The algorithm will recognize two parallel planes as the same patch, which can lead to errors in roof patch recognition. Meanwhile, because the method removes the building facade points, when the normal vector of the points is calculated, if the number of the points is considered to be more, the calculation of the normal vector is wrong.

in combination with this information, we propose a region growing algorithm in combination with the pitch of the points to segment the roof patch. Firstly, setting a threshold value of a point distance, clustering a point set by iteratively calculating the point distance of a nearest neighbor point, and counting the height mean value of each type of point. To prevent the effects of holes, we merge clusters of elevation approximations. Then, the different patches are identified by taking each type of point as input according to the principle of a classical region growing algorithm. Iterations are repeated until all points are traversed.

2.3 repair of roof shingles

Due to calculation errors and camera shooting angles, the image point cloud often has the situations that parts of points are inconsistent with the real world, even parts of points are missing, and the like, for example, a transition region between a roof plane and a vertical plane is an irregular curved surface, holes exist in a roof surface patch, and the like. Meanwhile, a roof of a building is generally provided with some accessories such as a water tank, a ventilation facility, a solar water heater, and the like. After facade removal, there will still be some residue. These conditions will cause the boundaries of the divided roof tiles to be irregular or even incomplete. Through observing a large number of building roofs, the building roofs are usually regular face sheets consisting of rectangles, triangles and sectors, points in the face sheets are distributed irregularly, but boundary points have strict shape constraints. Based on the point, the invention provides a new patch repairing method based on PCA, which can not only regularize patches, but also maintain the topological relation among the patches.

The algorithm mainly comprises two steps: the consolidation of roof shingles and appurtenances, and the repair of shingles. The algorithm principle is as follows:

1) All the dot flag bits are set to 0. Finding the patch P with the largest number of points, as shown in equation (9)_maxThe flag bit is changed to 1.

P_max＝{P_i|max{P_i)，i∈n} (9)

In the formula: p_iRepresenting the number of points in one panel and n representing the number of panels of the roof.

2) Calculating all patches and P_maxIs measured. In order to accurately find out the accessory point set, the identification position is 0 and 0.5<d<2.0. These patches are of the same type as P_maxThe relevant adjunct, collectively referred to as P'. Changing the identification bit of P' to 1;

3) Using principal component analysis algorithm to separate P' and P_maxSwitching to the main plane. Let the z coordinates of all points in P' be P_maxIs measured. And finishing projection.

4) On the principal plane, based on the mean dot spacingAnd (6) point supplementing is carried out.

i. Finding the maximum X of coordinates in a set of points_min、X_max、X_minand X_max；

Passing pointMaking a line L perpendicular to the x-axis_iWill be reacted with L_iAt a distance ofthe inner points are all projected on the straight line. Calculate the neighboring dot spacing dis ifIs inserted between two pointsAnd (4) points. Wherein i is1,2,3, …, up to X_min+i*d>X_max。

ii, with (0, Y)_min) Repeating step ii until Y is reached_min+i*d>Y_max。

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

note: in order not to destroy the topological relation of the boundary of the patch, a threshold value T is set₁E.g. T₁＝4.0。

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

6) And (5) iteratively executing the steps 1) to 5) until the identification bits of the points are all 1.

through the processing of the algorithm, a regular and complete roof patch of the building can be obtained theoretically.

2.4 refinement of roof tiles

It has been found through experiments that the image point cloud usually contains building interior points due to the camera shooting angle. In order to avoid identifying the interior points of the building as roof patches, a filtering method based on a two-dimensional grid is provided.

the principle is as follows: first, find the maximum X on the X, y axes_min、X_max、Y_min、Y_max(ii) a . And calculating the row-column index number of the grid to which each point belongs according to a formula (10). Counting the number of points in each grid, and sequencing the points according to the elevation values. If the difference between the elevations of neighboring points within the grid is greater than a threshold value Tz, such as Tz 2.5, the lower elevation portion of the points is discarded.

In the formula: row (Row)_icolumn index, representing point i_idenotes its column index, width denotes the width of the grid, x_iRepresenting the x-coordinate, y, of point i_iRepresenting its y coordinate.

Due to the adoption of the technical scheme, compared with the roof extraction technology of buildings based on airborne LiDAR, the roof extraction method has the following obvious advantages:

1) The image point cloud acquisition technology requires cheap and simple equipment and can obtain accurate and vivid models. The inclined imaging system mainly comprises a high-precision POS system and a plurality of cameras, and is low in manufacturing cost. On the contrary, the cost of the airborne laser radar equipment is high because one laser scanner needs hundreds of thousands and an inertial navigation system and other equipment are added. Meanwhile, the three-dimensional reconstruction technology based on the image can generate an accurate and vivid dense point cloud model, and the requirement of extracting the roof contour line of the building is met.

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

3) The full-automatic building roof contour line extraction algorithm can automatically acquire the roof contour line through simple parameter setting, thereby solving the problems of long time consumption and large manpower and material resource investment in traditional house cadastre measurement.

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 example 1;

FIG. 3 is a schematic view of the building as a unit in embodiment 1;

FIG. 4 is a schematic diagram of the point cloud of the roof patch in example 1;

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

Fig. 6 is a schematic diagram comparing the RANSAC segmentation algorithm and the method of the present invention.

Detailed Description

Reference will now be made in detail to embodiments of the present invention, examples of which are illustrated in the accompanying drawings, wherein like or similar reference numerals refer to the same or similar elements or elements having the same or similar function throughout. The embodiments described below with reference to the drawings are illustrative and intended to be illustrative of the invention and are not to be construed as limiting the invention.

In the description of the present invention, it is to be understood that the terms "center", "longitudinal", "lateral", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "clockwise", "counterclockwise", and the like, indicate orientations and positional relationships based on those shown in the drawings, and are used only for convenience of description and simplicity of description, but do not indicate or imply that the referenced devices or elements must have a particular orientation, be constructed in a particular orientation, and be operated, and thus are not to be considered as limiting the present invention.

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

In the present invention, unless otherwise expressly stated or limited, the terms "mounted," "connected," "secured," and the like are to be construed broadly and can, for example, be fixedly connected, detachably connected, or integrally connected; can be mechanically or electrically connected; they may be connected directly or indirectly through intervening media, or they may be interconnected between two elements. The specific meanings of the above terms in the present invention can be understood by those skilled in the art according to specific situations.

In the present invention, unless otherwise expressly stated or limited, "above" or "below" a first feature means that the first and second features are in direct contact, and that the first and second features are not in direct contact but are in contact with each other via another feature therebetween. Also, the first feature being "on," "above" and "over" the second feature includes the first feature being directly on and obliquely above the second feature, or merely indicating that the first feature is at a higher level than the second feature. A first feature being "under," "below," and "beneath" a second feature includes the first feature being directly under and obliquely below the second feature, or simply meaning that the first feature is at a lesser elevation than the second feature.

The technical solution of the present 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 the oblique images acquired by the unmanned aerial vehicle after the oblique images are subjected to a motion recovery structure algorithm is selected for experiment. The original data is collected in a certain cell of China, and a total of 21 txt files are different in size, each file comprises about 38593 ten thousand points, and the size of each file is about 17.6G. The types of buildings in the scene are complex, including 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 number of points in the original image point cloud is large and the point density is large, the data needs to be thinned. Meanwhile, a large number of separated block noise points exist under the ground in the scene through observation. The sparse and clustering algorithm is executed according to file iteration, after blocky noise is removed by the clustering algorithm, the sparse proportion is executed according to the ratio of 1:50, and finally the sparse and clustering algorithm is combined into a file, wherein the size of the file is 359M, and the point number of the file is about 772 ten thousand. As shown in fig. 2.

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

And a, separating ground points from 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:

7) calculating the distance between the particle and the point cloud;

for specific implementation, see the related documents Wuming Z, Jianbo Q, Peng W, et al, an Easy-to-Use air absorber LiDAR Data filtration Method Based on the same mechanism [ J ]. removal Sensing,2016,8(6): 501-.

Step b, utilizing the split G_rAnd identifying green or near-green color points by using the vegetation coefficient index VI, and optimizing the identification point set by using a change detection algorithm based on a normal vector so as to eliminate the vegetation points.

In the embodiment, in the step b, because the vegetation is generally green or near green, the DN (Digital Number) value of the green light band of the vegetation is generally higher than the DN values of the red and blue light bands. In combination with these two points, the present invention proposes the use of the split G_rAnd identifying vegetation according to the vegetation coefficient index VI. The definition of the split and vegetation coefficient is shown as follows:

By setting the split G_rAnd a threshold value of the vegetation coefficient index VI, such as G_r>0.35，VI>0, a green or near-green set of points can be effectively identified.

However, buildings of special colors are often found in the real world, such as kindergartens with greenish-looking walls in the examples. In order to prevent the vegetation points from being mistakenly removed during filtering, a change detection algorithm based on a normal vector is adopted to further screen an identification point set and determine the vegetation points.

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

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

the maximum change in the normal vector can be quantified using eigenvalues of the covariance matrix based on normal vectors as mentioned in M.Pauly, Point principles for Interactive Modeling and Processing of 3D geometry.Konstanz, Germany: Hartung-Gorre,2003.

in particular, neighborhoods are constructedthe covariance matrix of the normal vectors of all points in the block is shown in equation (3). Same as C_pCalculating eigenvalues of the covariance matrix, sorting the eigenvalues,The maximum change of the normal vector of the point p on the gaussian sphere can be quantitatively represented.

By setting a threshold value Tⁿand non-vegetation points concentrated by near-green color points can be effectively identified.

And c, performing connected domain analysis by using an Euclidean clustering algorithm to obtain a preliminary building monomer result, and further optimizing by using a constraint term based on an elevation threshold value.

Example step c, the point cloud obtained by step b typically comprises buildings (including fences), vehicles and power implementations such as wires, poles, etc. Connected domain analysis is carried out through the Euclidean clustering algorithm, most non-building points can be proposed by setting a point distance threshold value and a minimum point number threshold value, and preliminary building clustering is realized.

Specifically, the european clustering algorithm is a clustering algorithm based on a point distance, and the algorithm implementation steps are as follows:

1) Setting a threshold T of dot spacing, all dots being set as unaccessed dots;

4) repeat execution 2) -3) until all points are marked as access points.

after the euclidean clustering algorithm is executed, some non-building points are still identified as buildings in the embodiment, such as residual ground, plane-like green vegetation, walls connecting the buildings, and the like, and the presence of these misclassifications affects the accuracy of roof patch extraction. To remove these misclassifications, a constraint term based on an elevation threshold needs to be set here.

Specifically, the minimum elevation value Zmin of each type of point set is found, an elevation threshold value is set by combining building roof common knowledge, for example, z is greater than Zmin +3, and most non-building point sets can be filtered. After that, the euclidean clustering algorithm is performed again. Through these steps, building singleness with high accuracy can be realized. As shown in fig. 3 as a building unit.

And taking the point cloud of the single building as a processing unit, and executing subsequent processing steps.

And d, removing the vertical face of the building by using the point characteristics f based on the normal vector to obtain the top point of the building house.

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

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

f＝1-|n_p·e_z| (4)

in the embodiment, the point features f can only eliminate the wall points of the building, and can not be effectively removed for non-wall points, such as structures like a canopy and the like. And removing the data by adopting an Euclidean clustering algorithm.

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

In the embodiment step e, the roof of the building is generally composed of various planes or plane-like planes, which provides a thought for people. The roof patch 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, i.e. by setting thresholds for curvature variation and normal angle difference, merging points close in smoothness, and thus dividing a patch. Thus, each time the algorithm outputs, it is considered to be a co-planar point. However, this approach has some disadvantages here. The algorithm will identify two parallel planes as the same patch, which can lead to errors in roof patch identification. Meanwhile, because the method removes the building facade points, when the normal vector of the points is calculated, if the number of the points is considered to be more, the calculation of the normal vector is wrong.

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

in step f, the roof panel obtained in step e usually has irregular or incomplete boundaries, and even holes. By executing the PCA-based patch repairing method, holes can be repaired, patches can be flattened, and the topologic relationship among the patches can be kept.

Specifically, the algorithm mainly comprises two steps: the consolidation of roof shingles and appurtenances, and the repair of shingles. The algorithm principle is as follows:

P_max＝{P_i|max{P_i)，_i∈n) (9)

2) calculating all patches and P_maxis measured. In order to accurately find out the accessory point set, the identification position is 0 and 0.5<d<2.0. These patches are of the same type as P_maxThe relevant adjunct, collectively referred to as P. Changing the identification bit of P to 1;

3) Using principal component analysis algorithmsP and P are_maxSwitching to the main plane. Let the z coordinates of all points in P be P_maxIs measured. And finishing projection.

i. Finding the maximum X of coordinates in a set of points_min、X_max、Y_minAnd Y_max；

Passing through a pointMaking a line L perpendicular to the x-axis_iwill be reacted with L_iat a distance ofThe inner points are all projected on the straight line. Calculate the neighboring dot spacing dis ifIs inserted between two pointsand (4) points. Wherein i is 1,2,3, …, up to X_min+i*d>X_max。

iii, with (0, Y)_min) Repeating step ii until Y is reached_min+i*d＞Y_max。

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

5) and (5) restoring the P from the main plane to the world coordinate system to obtain the roof patch P.

And (5) iteratively executing the steps 1) to 5) until the identification bits of the points are all 1.

And g, utilizing the two-dimensional grid to realize the duplicate removal treatment of the roof surface patch for the output of the step e.

In the embodiment, in step g, the output patches obtained in step f are partially building internal patches due to the camera shooting angle and the building window. For this purpose, a filtering method based on a two-dimensional grid is adopted to eliminate the part of the point set.

The principle is as follows: first, find the maximum X on the X, y axes_min、X_max、Y_min、Y_max. And calculating the row-column index number of the grid to which each point belongs according to the formula (10). Counting the number of points in each grid, and sequencing the points according to the elevation values. If the difference between the elevations of neighboring points in the grid is greater than a threshold value Tz, for example, Tz is 2.5, the lower elevation part of the points is discarded.

the result of dividing the roof of the building into patches in this embodiment is shown in fig. 4. Experiments show that the partitioning algorithm can effectively partition a roof surface structure formed by combining rectangles in a data set. The results of the segmentation were quantified by experimental data and are counted below.

TABLE 1 roof segmentation experiment structure statistics

Wherein BN represents the number of the building; SP, representing the data of roof surface slices obtained after the roof of the building is divided; RP, which represents the number of roof tiles inherent to the building in the data set; complete, namely the ratio of the number of divided roof patches to the number of inherent roof patches;

The experimental result shows that the roof surface patch segmentation method provided by the invention can obtain better segmentation results in most buildings, the algorithm stability is good, and the average segmentation result accuracy is 97.37%. Comparing the experimental result with the original data, finding that the reason that the f and g groups of experimental results are under-segmented is that the two roof surface patches are positioned on the same horizontal plane and adjacent to each other, so that the roof surface patches cannot be segmented; the m and n groups of experimental data are similar in structure, the reason for the insufficient segmentation is that the separation distance between two roof surface patches is too small, so that the surface patches are combined, and the segmentation is insufficient due to the small area and too few points of the surface patches; the reason for over-segmentation of the h group is that a large hole is formed in the dough sheet and the number of attachments nearby is too large, so that after the building facade is removed, the area of the hole is increased, and over-segmentation is caused. These experimental results are sufficient to demonstrate that the algorithm of the present invention has good adaptability to building roof structures formed by combinations of convex polygons.

When the system is specifically implemented, the corresponding system is realized in a modularized mode. The system specifically comprises the following modules:

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

The module 2 is used for obtaining initial vegetation points by utilizing the green signal ratio and the vegetation coefficients and detecting and further restraining the vegetation points based on the change of the normal vector;

The module 3 is used for executing an Euclidean clustering algorithm and elevation threshold value constraint so as to realize building singleness;

The module 4 is used for realizing a building facade elimination algorithm and realizing roof surface patch segmentation to obtain an initial roof surface patch;

and the module 5 is used for realizing a patch repairing algorithm based on PCA and patch duplication removal based on the two-dimensional grid to obtain a flat roof patch.

the specific implementation of each module can participate in the corresponding step, and the detailed description is omitted here.

The method provides a framework flow and relates to linear combination of a plurality of algorithms. As shown in tables 2 and 3, the number of parameters on the whole frame surface is up to 18, and actually, the number of parameters to be adjusted is only 4, which are the threshold Tf of the point feature f, the flatness threshold SmoothnessThreshold and curvature threshold CurvatureThreshold in the region growing algorithm, and the neighborhood point number or radius, respectively. Each parameter is analyzed in detail below.

TABLE 2 parameters of building singulation technique

TABLE 3 parameters in roof tile segmentation technique

In table 2, the vegetation filtering section filters green plants based on the color information of the dots, so Gr and Vi are constants; threshold value of variation T for point normalⁿIt was found in a number of experiments that after normalization of the mean-based values, TⁿThe ideal effect is always obtained when the value is 0, and thus it is considered that the value is constant. And the cloth filtering part can meet the requirement by using default parameters. Because by setting the elevation threshold T in the singulation technique_zThe ground filtering results can be optimized. In singulation, the height of the building is always greater than 3m in combination with the common sense, so that the elevation threshold T is_z2.5 is a constant; the minimum point number minP and the dot pitch threshold dis are usually constant values, and unless the building specification or the dot density is largely changed, minP is 10000 and dis is 0.5, respectively.

In table 3, the threshold T1 for the facade removal section, point feature f, is typically equal to approximately 1, where appropriate values need to be set by testing of the sample building. For the Euclidean clustering algorithm, both the singularization and the patch segmentation are involved, the threshold value of the point distance is 0.5, and the minimum point number minP is changed here, usually about 500. The region growing algorithm in the patch segmentation needs to properly adjust two parameters, namely SmoothnessThreshold and Curvaturethreshold according to the quality of the image point cloud. A patch repairing part, Td is mainly used for preventing the topological structure of a building from being damaged, and generally has a minimum value according to the common knowledge of building construction, namely Td is 3; dot pitch threshold T₁usually a constant value, T₁＝0.05. In the patch optimization section, the grid width is typically set to 3 times the dot spacing.

It is noted that the algorithm involves the estimation of the normal many times, and the method adopts a neighborhood-based principal component analysis method. Because the coordinates of the image point cloud points are shifted compared with the real situation, the neighborhood setting must be larger, usually the number of nearest neighbor points is set to 500 or the radius of the neighborhood is about 0.5 m, and the adjustment needs to be made according to the quality of the image point cloud.

The method is a combined application of a series of algorithms, and a plurality of sub-algorithms are involved in the method. In the prior art, the references only relate to partial principle overlapping, such as a region growing algorithm, a Euclidean clustering algorithm, two-dimensional grid formation, a principal component analysis method and the like.

The following is a detailed analysis of each comparison file to find out the prior art with real comparison significance.

CN107220987A records a rapid detection method for the roof edge of a building based on principal component analysis, which aims to utilize a principal component analysis algorithm to assist in detecting the edge of a roof patch, and the invention utilizes the principal component analysis algorithm to repair the roof patch. The key points of the invention are obviously different.

CN108010092B records an urban high-density area solar energy utilization potential evaluation method based on low-altitude photogrammetric data modeling. In this document, in the building singleton part, it attempts to separate the building point cloud from the color information in the image point cloud by means of a volvox plug-in a parametric modeling tool, gradsphopper. It has the disadvantage that it does not take into account the fact that the colour of a building is similar to that of the ground or vegetation. The key points of the invention are obviously different.

CN107944384A (covariance matrix-area growth) discloses a building roof tile segmentation method based on three-dimensional Voronoi diagram, and the processing object is an airborne LiDAR point cloud. According to the method, the three-dimensional Voronoi graph is used for replacing kdtree to establish the point neighborhood, parameter setting during neighborhood setting is avoided, and the problem that the neighborhood is not easy to control during point cloud data space topological relation construction is effectively solved. However, because the processing object is an airborne LiDAR point cloud, the vegetation filtering algorithm adopted by the processing object is completely different from the method and has obvious difference from the core key point of the method.

CN107230251A discloses a technique for creating a 3D city model from oblique imaging data and lidar data. This document relates to a method of matching fusion of oblique imaging data and lidar data to create a 3D city model and the design of a hybrid 3D imaging device. The invention describes an algorithm for processing colored three-dimensional point clouds, which have no contrast. The key points of the invention are obviously different.

CN108090957A discloses a BIM-based method for mapping terrain. The document introduces a method for establishing a three-dimensional point cloud model mainly by human-computer interaction, and the invention relates to a full-automatic building roof surface patch identification algorithm, which is obviously different from the core key point of the invention.

CN106846494A (grid method _ GPU) discloses an automatic simplex algorithm for oblique photography three-dimensional building models. The algorithm transmits vertex information of the building contour line into a GPU shader, and highlights the building in the contour line by marking the building in the contour line, so that the building is subjected to building singleization operation. The document adopts the GPU technology, and the technical route of the GPU technology is greatly different from that of the invention. The key points of the invention are obviously different.

CN108074232A (grid-voxel) is the same as CN108109139A, disclosing an onboard LiDAR building detection method based on voxel segmentation. The method comprises the steps of regularizing original airborne LiDAR point cloud data into a gray 3D volume element data set, and dividing and marking the gray 3D volume element data into a plurality of 3D connected regions based on connectivity and radiation characteristic similarity criteria; voxel segmentation-based on-board LiDAR building detection is then completed based on the characteristics of the building rooftops and facades. The method adopts a voxel technology and is obviously different from the core key point of the invention.

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

CN102411778B discloses an automatic matching method for airborne laser point cloud and aerial image. The method directly extracts a building contour line from an airborne LiDAR point cloud, obtains the building angle characteristics of registration primitives according to the building contour line, and then automatically matches the homonymous angle characteristics between the point cloud and an image under the assistance of an aerial image approximate exterior orientation element; and then, adopting a beam method block adjustment and a loop iteration strategy to realize the integral optimal registration of the aerial image and the point cloud data. The key points of the invention are obviously different.

CN101726255B discloses a method for extracting a building of interest from three-dimensional laser point cloud data. The method attempts to directly separate the building from the non-ground point by using the Euclidean clustering algorithm, and theoretically, the method has certain limitation. Because non-building points such as vegetation and large automobiles are always on the ground in a real scene, the robustness of the simple execution of the Euclidean clustering algorithm is not high. Meanwhile, the purpose of this document is to search for buildings by setting building edge features. Therefore, the method is obviously different from the core point of the invention.

CN107545602A (spatial topological relation-building modeling) discloses a building modeling method under the spatial topological relation constraint based on LiDAR point cloud. The document focuses on spatial topological relation processing of building rooftop geometric elements and has little detailed description on building singulation and rooftop segmentation. The key points of the invention are obviously different.

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 patch, the roof patch of the building is simple in segmentation part, the roof patch is obtained by directly adopting a region growing algorithm, and then a plane equation of a plane where the roof patch is located is obtained by fitting with a least square method. The key points of the invention are obviously different.

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

The method for extracting trees based on region growing and gradient segmentation is provided for extracting trees from urban LiDAR point clouds, and has no contrast with the method.

CN106970375A discloses a method for automatically extracting building information in an airborne laser radar. The method essentially detects all building point cloud planes by using a plane fitting method, and then removes part of misclassified point clouds by using three-dimensional morphological corrosion operation, thereby obtaining accurate building point clouds. The method is suitable for airborne LiDAR point cloud, however, the image point cloud is different from the laser point cloud, the density of the point cloud is high, and the method also comprises complete building vertical plane points while the roof point cloud is contained. Therefore, if no other measures are taken, the plane fitting method is directly carried out, and the building point cloud cannot be detected through a small amount of iteration. Therefore, the method is obviously different from the core point of the invention.

a building roof point 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 in the traditional RANSAC algorithm to a certain extent. The patch segmentation algorithm of the invention adopts a region growing algorithm, and has obvious difference with the core key point of the invention.

Journal paper 'improved RANSAC point cloud segmentation algorithm considering building roof structure' utilizes a gradient and height difference triangular region growing method to decompose different structural levels of a complex building, and then provides a RANSAC algorithm with a floating congregation threshold value to detect a building roof plane, so that the method has certain applicability. The patch segmentation algorithm of the invention adopts a region growing algorithm, and has obvious difference with the core key point of the invention.

A research graduate thesis 'airborne LiDAR point cloud data filtering and building point group segmentation research' improves a method for extracting point clouds contained in a building roof surface patch by random Hough transform. The assumption of the whole algorithm is that the building model is composed of a plurality of planes, and the assumption has certain limitation because in a real scene, part of buildings are composed of curved surfaces. The key points of the invention are obviously different.

CN109242855A discloses a roof segmentation method, system and device based on multi-resolution three-dimensional statistical information. The method comprises the steps of extracting three-dimensional point cloud feature statistical information with different resolutions from an image point cloud, and then performing semantic classification on a three-dimensional point cloud scene by directly utilizing global energy optimization to obtain the point cloud of a building. The whole technical route is different from the method adopted by the invention.

CN108171720A discloses a method for detecting the object boundary of oblique photography model based on geometric statistical information. The method includes the steps of calculating a two-dimensional minimum outsourcing rectangle of the building, enabling the minimum outsourcing rectangle to be intersected with a model, and selecting point cloud data belonging to an area where a current object is located to obtain a single building. The method has complex process and large calculation amount. The method is different from the technical route adopted by the invention.

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

CN109461207A discloses a method and device for building singleization by point cloud data. Before the classification processing is carried out on the voxels, the method carries out preprocessing on the three-dimensional points in the voxels, and divides the three-dimensional points 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. The method is greatly different from the technical route adopted by the invention.

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

CN109859315A discloses a method for separating earth surface vegetation in three-dimensional images. The method is a software operation method and has no contrast with the method.

CN109870106A discloses a building volume measuring method based on unmanned aerial vehicle pictures. According to the method, building image point cloud is obtained through the dense matching of unmanned aerial vehicle pictures, and then the volume of a building is calculated for a model after Delaunay triangularization by using an integral method. The key points of the invention are obviously different.

CN105844629B (RGA-DP-RANSAC) discloses a large-scene city building facade point cloud automatic segmentation method. The method for extracting the point cloud data of the roof 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 triangulation network encryption, then filtering out feature points with the height difference smaller than 2.0m by taking the elevation of the ground points as a reference, and then segmenting the roof point cloud of the building by adopting an RANSAC patch detection algorithm. The method is relatively simple and is not suitable for image point clouds with a large amount of elevation information. The key points of the invention are obviously different.

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 point cloud echo times are not available in the image point cloud, so that the method is not applicable. Secondly, it adopts traditional region growing algorithm to identify the building point cloud, which has certain limitation. And 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 point cloud space distribution characteristics, optimization is performed by utilizing a patch normal vector included angle constraint and a patch distance constraint, and similar patches are combined. The key points of the invention are obviously different.

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

CN107644452A (onboard LiDAR-roof patch segmentation) discloses an onboard LiDAR point cloud roof patch segmentation method and system. The method is different from the invention in that:

1) The method includes the steps that a point cloud neighborhood system is built by taking a 3D Voronoi diagram as airborne LiDAR point cloud data, and the 3D Voronoi neighborhood-based adjacency relation is built for any point in the airborne LiDAR point cloud data. The invention employs kdtree.

2) This document uses multiple echo information from a LiDAR point cloud to distinguish between building points and vegetation points. The invention processes data into image point cloud without echo information, and adopts other methods.

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

1) And filtering vegetation in the ground points by adopting color invariants, and not considering the condition that buildings have similar colors. After the color filtering is adopted, the change detection algorithm based on the normal vector is added, so that the color filtering is ensured not to cause the loss of building points.

2) The building point cloud is obtained by using a region growing method based on plane fitting, and is essentially based on the assumption that the building is composed of planes. After vegetation and the ground are filtered, the method directly adopts the Euclidean clustering algorithm to obtain the building by setting a point interval threshold and the minimum clustering point number. It is not based on this assumption, ensuring that buildings of arbitrary structure can be identified.

3) And a facade dividing section which separates the facades of the building by using the direction of the normal vector of the patch. The present invention uses normal vectors of points.

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

1) Grasping a method that the density of the facade points of the building in the airborne LiDAR point cloud is smaller than the roof top point to put forward the wall points;

2) The roof layer resampling part uses a support vector machine algorithm, and the calculation amount is increased.

3) the method comprises the following steps of selecting a seed area, growing a roof surface patch, flattening and optimizing the surface patch, obtaining initial plane parameters by using a point curvature threshold preliminarily, and then realizing the growth of the roof surface patch by using a distance threshold and a distance standard deviation.

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

In terms of building identification, an elimination method is employed. Namely, the land features are classified into 3 types: building, ground and vegetation, and removing ground points and vegetation points, wherein the rest are building points. Many of the methods in the references are based on the assumption that a building is composed of planar patches, which has certain limitations and cannot extract a building composed of non-planar patches. The present invention does not have this problem.

In the aspect of roof surface patch segmentation, after the original roof surface patch point cloud is extracted, a surface patch repairing method based on PCA is adopted, so that the problem that holes exist in the original image point cloud is solved, and the surface patch is flattened. In the reference, there is no reference to the processing of holes in the image point cloud.

In particular, the method comprises the following 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 the principle and the like, and meanwhile, the method has strong scene applicability, not only can be used for processing flat ground, but also is suitable for the ground with a certain gradient. While the majority 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 variance of the normal vector in the green points is measured by eigenvalues in the covariance matrix to determine which points belong to the building points. Under the dual action of the front and back steps, the green vegetation can be removed, and meanwhile, the building points can not be deleted by mistake. In the reference, there are methods for filtering vegetation by using color information, but they do not consider the condition that the color of a building is similar to the vegetation, so the algorithm has a certain limitation.

3) In a roof patch repairing algorithm, a patch repairing method based on PCA projection conversion is provided, which can process attachments of a roof, flatten patches, repair holes and keep topological relation among patches. The reference does not consider the problem of repairing holes in the wafer.

4) In the refinement part of the roof surface patch, a method for identifying overlapped surface patches based on a two-dimensional grid is adopted. In airborne LiDAR point cloud or unmanned aerial vehicle slope image matching point cloud, because the inclination angle problem of laser radar or camera, often contain some building interior points or even the inside dough sheet point sets of building in the point cloud, here through carrying out the elevation to graticule mesh point set and gathering the class, the mode of rejecting the point to concentrate the height lower person, accomplish the deletion of building interior point to the accuracy of building roof dough sheet discernment has been improved.

the following demonstrates a comparison of the method of the invention with a prior art example.

1. Vegetation filtering algorithm aspect:

the plant filtering method based on color invariant as mentioned in CN104484668B and the method of the present invention were compared, and the data set was tested with the data set in the specific examples.

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 for the color invariant be Tg, and the color invariant formula defined by the green and blue color 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_g(x, y, z) represents the color invariant at the (x, y, z) point. When psi_g＜T_gWhen the point is a vegetation point, the point is represented as a vegetation point; otherwise, it is a non-vegetation point.

TABLE 3 Vegetation Filter Algorithm comparison

Method of producing a composite material	Pts
		The method of the invention	216335
theory of color constancy	225893

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

2. And (3) roof surface patch segmentation algorithm aspect:

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

TABLE 4 roof partition experimental structure statistics

In order to better compare the two algorithms, three buildings of g, m and n are selected for individual analysis. As shown in fig. 6, the results are extracted for the roof tiles of three buildings, g, m, and n. Wherein (a-c) is a RANSAC method and (d-f) is the method of the present invention. Results show that the RANSAC point cloud segmentation algorithm is sensitive to the number of patch points and cannot effectively extract small patches. Therefore, the algorithm of the invention has better stability and can better keep the detailed structure of the roof surface.

In the description herein, references to the description of the term "one embodiment," "some embodiments," "an example," "a specific example," or "some examples," etc., mean 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, the schematic representations of the terms used above do not necessarily refer to the same embodiment or example. 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 is understood that the above embodiments are exemplary and should not be construed as limiting the present invention, and that variations, modifications, substitutions and alterations can be made in the above embodiments by those of ordinary skill in the art without departing from the principle and spirit of the present invention. Any modification, equivalent replacement, improvement and the like made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. a method for segmenting a roof surface patch of a building from a large-scale image dense matching point cloud is characterized by comprising the following steps:

3) Removing the building facade;

4) And (4) dividing, repairing and refining the roof surface patch.

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

3. The method for segmenting a roof patch of a building from a large-scale image dense matching point cloud according to claim 1, wherein the ground filtering is implemented as follows:

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

4. The method of claim 3, wherein the ground filtering is implemented by a post-processing step comprising: searching four adjacent particles of each movable particle, and if finding that the immovable particles exist, comparing the elevations of the movable particles and the immovable particles to the stress light spot; if the minimum height difference is less than the threshold, the movable particle continues to move until set to immovable; this step is performed iteratively until all movable particles are traversed.

5. The method of claim 1, wherein the vegetation filtering utilizes the split green ratio G_rAnd identifying green or near-green color points by using the vegetation coefficient index VI, and optimizing the identification point set by using a change detection algorithm based on a normal vector so as to eliminate the vegetation points.

6. The method of claim 5, wherein the Luxin ratio is the ratio of DN values of green bands to the sum of DN values of RGB; the definition of the split and vegetation coefficient is shown as follows:

7. The method as claimed in claim 1, wherein the building is singulated, connected domain analysis is performed by using Euclidean clustering algorithm to obtain preliminary building singulation result, and then further optimization is performed by using a constraint term based on elevation threshold.

8. The method of claim 1, wherein the removing of the building facade is performed by removing the building facade using a point feature f based on a normal vector to obtain a building roof point cloud; the definition of the point feature f is shown as follows:

f＝1-|n_p·e_z|

9. The method for segmenting a building roof patch from a large-scale image dense matching point cloud as claimed in claim 1, wherein the segmentation of the roof patch is to perform the segmentation of the roof patch by using a region growing algorithm based on a point spacing for a single building roof point cloud; the repairing of the roof surface patch refers to that the hole is repaired by utilizing a patch repairing algorithm based on PCA (principal component analysis) on the output of the divided roof surface patch, and the patch is flattened; the refinement of the roof surface patch refers to that the two-dimensional grid is utilized to realize the duplicate removal processing of the roof surface patch for the output of the divided roof surface patch.

10. a system for segmenting a building roof tile from a large-scale image dense matching point cloud using the method of any one of claims 1 to 9, comprising: