Window detection method based on hole constraint and hierarchical positioning

Abstract

The application discloses a window detection method based on hole constraint and hierarchical positioning, first extracts the building facade in a complex outdoor scene; then with the help of the hole idea, the main wall surface extracted from the building facade is subjected to hierarchical judgment to obtain accurate window boundary information, and the window boundary point is individualized through clustering; then for each single window boundary point, the feature points thereof are extracted, and the boundary points are regularized in the mode of equidistance interpolation points; finally, a window line frame model is constructed to obtain a window detection result. The application solves the problem of inaccurate window detection caused by irregular distribution and non-uniform shape of windows in a complex facade, and can realize window detection of complex buildings with different facade structures, and improves the detection accuracy.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of computer vision and image processing method, and relates to a window detection method based on hole constraint and hierarchical positioning. BACKGROUND

[0002] With the increasing application in the fields of city planning, digital city construction, virtual tourism and the like, the demand for three-dimensional building reconstruction is also increasing. Among them, the detailed building facade structure is an important part of three-dimensional building reconstruction. As the main component of the building facade, the window is the key to the construction of the fine facade model.

[0003] Due to the reason that the glass material of the window in most building facades cannot reflect the laser beam, and the complexity of the outdoor scene, such as the shielding of the elements such as vegetation, the accuracy and integrity of the window detection are seriously affected. In addition, most of the current researchers mainly focus on the modern building facade in urban areas, and detect the rectangular window with repeated structure and symmetry characteristics. There are few related researches on the window detection of irregular building facades and low-rise residential buildings in rural areas (windows with various shapes, sizes and asymmetric positions on the facade). SUMMARY

[0004] The purpose of the application is to provide a window detection method based on hole constraint and hierarchical positioning, which solves the problem of window detection of different facade structures of complex facades.

[0005] The technical scheme adopted by the application is that the window detection method based on hole constraint and hierarchical positioning is implemented according to the following steps:

[0006] Step 1: filtering and rotating the point cloud original data to rotate the scene to be parallel to the ground and the horizontal plane;

[0007] Step 2: on the basis of step 1, removing the ground points in the scene, clustering the non-ground points, and extracting the building points by analyzing the characteristics of the non-ground points;

[0008] Step 3: extracting the required building facade;

[0009] Step 4: extracting the main wall surface in the extracted building facade;

[0010] Step 5: on the main wall surface, extracting the window boundary points by using the hierarchical method;

[0011] Step 6: performing window boundary optimization operation.

[0012] The application also has the characteristics that,

[0013] Step 1 is specifically:

[0014] Step 1.1: mixed filtering operation on original point cloud data

[0015] Step 1.1.1, divide the original point cloud data into several cubes with edge length l, calculate the center of gravity of the cube and use the center of gravity point to represent all points in the cube, and replace the center of gravity point with the nearest point in the original point cloud data;

[0016] Step 1.1.2, continue to remove noise points in the form of radius filtering: judge whether the number of nearest neighbor points in the specified radius sphere is less than the specified threshold as the basis for whether the point is a discrete point, if the number of nearest neighbor points in the specified radius sphere is less than the specified threshold, it is determined that the point is a discrete noise point and is deleted; otherwise, it is retained;

[0017] Step 1.2: realize scene rotation

[0018] Step 1.2.1, according to the feature that the ground point normal vector is parallel to the y axis, the points with normal vector and y axis included angle less than are preliminarily extracted, and plane fitting and the normal vector n1 of the plane are obtained;

[0019] Step 1.2.2, give the normal vector n2 of the target plane;

[0020] Step 1.2.3, calculate the rotation angle and rotation axis: calculate the included angle θ of n1 and n2 according to formula (1), that is, the rotation angle; the rotation axis is the vector perpendicular to the normal vector n1 of the main plane and the normal vector n2 of the target plane,

[0021] θ = arccos(n1, n2) (1);

[0022] Step 1.2.4, solve the rotation matrix R according to the rotation axis and the included angle θ, and transform the original point cloud scene through the rotation matrix R.

[0023] In step 1.2.1, is 10°.

[0024] The step 2 is specifically as follows:

[0025] Step 2.1, reverse the point cloud data obtained in step 1, generate an initial cloth grid according to the set grid size, and determine the number of grid particles, wherein the starting position of the cloth is usually above the highest point of the reversed point cloud;

[0026] Step 2.2, project the point cloud data and the grid particles to the same horizontal plane, search for the nearest neighbor point of each grid particle, and record the elevation value IHV before projection;

[0027] Step 2.3: For each grid particle, if it can move, calculate its displacement caused by the external force and compare the position of the particle after displacement with the IHV of its corresponding nearest neighbor. If the height of the particle is less than or equal to the elevation value, the particle is considered to be an immovable point and its height is set to IHV.

[0028] Step 2.4, calculate the displacement distance of grid particles caused by internal factors;

[0029] Step 2.5: Repeat steps 2.3 and 2.4 until the maximum number of iterations is reached or the moving distance of all particles is less than the preset value, then stop the iteration;

[0030] Step 2.6, calculate the distance between the grid particle and its nearest neighbor: if it is less than the distance threshold dis g , it is considered to be a ground point, otherwise it is a feature point. In this way, the preliminary ground point set G can be obtained;

[0031] Step 2.7, for the set G, calculate its average height and the normal vector of each point;

[0032] Step 2.8: For each point in the set G, if the height of the point is higher than the average height, calculate the angle between its normal vector and the z-axis. If the angle is less than the set threshold α, skip it; otherwise, classify it as a non-ground point set and delete it from the ground point set.

[0033] Step 2.9, search the number of non-ground points around the ground seed point through K-nearest neighbor (KNN);

[0034] Repeat steps 2.8 and 2.9 until all seed points have been traversed.

[0035] The step 3 is specifically as follows:

[0036] Step 3.1: Extract buildings by fusing dimensionality features and roughness.

[0037] Step 3.2: Use the improved RANSAC method to further extract the building facade.

[0038] Step 3.1 is as follows:

[0039] Step 3.1.01, input the non-ground point set Q obtained in step 2 = {q1,q2,…,q m}, m is the total number of non-ground points, and a kd tree is created to find the nearest neighbor points N of the search point q through KNN q ; Next, calculate N q The covariance eigenvalues ​​λ1, λ2, λ3 of the inner point cloud (λ1≥λ2≥λ3); finally, the surface index λ is calculated by formula (2)P , divergent index λ S and rod index λ L , collectively referred to as dimension features, i.e. calculating the dimension feature of each point q j

[0040]

[0041] Step 3.1.02, label all points in set Q with label, set label, and initialize the label of the cluster segmentation block as segLab;

[0042] Step 3.1.03, traverse all points, if the label of the current seed point is equal to the set value, then find the nearest associated point q NN by searching for neighboring points, and store it in array V q , filter the associated points by dimension threshold, delete the associated points with different features from the seed point from V q , and thus obtain a new associated point array V n ;

[0043] Step 3.1.04, traverse all points in array V n , if there is a point whose label is not equal to the set value, find the smallest label value, compare it with segLab, and take the minimum value of the two as minSegLab; if there is no point whose label is not equal to the set value, directly assign segLab to minSegLab;

[0044] Step 3.1.05, traverse array V n again, if the label of the current point is greater than minSegLab, traverse all points again, if the label of a point is the same as that of the current point, reset the label value of the point to the value of minSegLab;

[0045] Step 3.1.06, add 1 to the value of segLab, repeat steps 3.1.3-3.1.5 until the value of label of all points is not equal to the set value, thereby completing the rough segmentation of the entire point cloud;

[0046] Step 3.1.07, calculate the average height and roughness R s of each segmentation block F i (i∈1,2,…,N) by formulas (3) and (4) respectively;

[0047]

[0048] ​

[0049] In formula (3) and (4): - average height of the segmentation block; Z j - height of each point in the segmentation block; n - number of points in the segmentation block; R s - roughness; N - number of segmentation blocks;

[0050] Step 3.1.08, regarding each segmentation block as a whole, adding a classification label to it, and taking the first segmentation block as a seed surface, starting from the seed surface to query;

[0051] Step 3.1.09, finding all adjacent segmentation blocks of the surface through spatial connectivity analysis. When the similarity threshold of the seed surface and the adjacent segmentation block is different, it can be considered that the segmentation block categories of the two are different, and merging cannot be performed, then the next adjacent segmentation block is selected; otherwise, the label of the adjacent segmentation block is reset to the label of the seed surface, and the next adjacent segmentation block is judged. The similarity threshold is set according to the roughness of the point cloud surface to judge the similarity of different segmentation blocks;

[0052] Step 3.1.10, judging whether all adjacent segmentation blocks have been traversed, if there are still adjacent segmentation blocks that have not been traversed, repeating step 3.1.9 until the merging of all adjacent segmentation blocks is completed; otherwise, the segmentation block that has not been accessed is taken as a seed surface;

[0053] Step 3.1.11, judging whether all segmentation blocks have been accessed, if there are still segmentation blocks that have not been accessed, repeating steps 3.1.9-3.1.10; otherwise, ending the merging;

[0054] Step 3.1.12, after the merging is completed, classifying the point cloud data according to the classification label;

[0055] Step 3.1.13, judging the number of points of each class cluster after merging, setting a maximum number threshold max and a minimum number threshold min, and retaining the class clusters within this range;

[0056] Step 3.1.14, combining different ground features, using elevation and geometric features to screen the retained class clusters, extracting the class clusters belonging to building points, and realizing accurate extraction of buildings.

[0057] Step 3.2 is specifically:

[0058] Step 3.2.1, in the building point cloud data obtained in step 3.1, randomly selecting multiple points to form a plane, and calculating the parameters of the plane model;

[0059] Step 3.2.2, calculating the distance of all points in the data set to the plane, if the distance is less than a distance threshold dW If the point is within the range, it is classified as an inlier, otherwise as an outlier;

[0060] Step 3.2.3, repeat steps 3.2.1-3.2.2 until the iteration number reaches the iteration threshold, then extract the plane with the largest number of inliers;

[0061] Step 3.2.4, count the total number of points in the plane, if it is less than the threshold S num , then delete it, otherwise project the extracted plane to the YOZ plane, so that all points form a vertical line segment, sort the projection points according to the y value from small to large, and use formula (5) to calculate the distance between adjacent two points;

[0062]

[0063] In formula (5):

[0064] P i —points in facade W j , j = 1, 2, …, s;

[0065] K—number of points in facade W j ;

[0066] Step 3.2.5, determine whether there is a vertical gap in the facade by determining which points have an association threshold greater than 3 times the median distance threshold Dis y , where each gap can divide the data into two different point clusters, if there is a gap, then the facade is divided again to optimize the extraction result, otherwise skip;

[0067] Step 3.2.6, project the extracted plane to the XOZ plane, so that all points form a horizontal line segment, sort the projection points according to the x value from small to large, and use formula (5) to calculate the distance between adjacent two points;

[0068] Step 3.2.7, determine whether there is a horizontal gap in the facade by determining which points have an association threshold greater than 3 times the median distance threshold Dis x , where each gap can divide the data into two different point clusters, if there is a gap, then the facade is divided again to optimize the extraction result, otherwise skip;

[0069] Step 3.2.8, repeat steps 3.2.1-3.2.7 until all potential facades are extracted from the building point cloud, then extract the required building facade from all potential facades according to the main wall surface characteristics, including height, perpendicular to XOZ plane characteristics, and add it to the set W, W = {W1, W2, … W s}, where W jThe facade is denoted as the jth facade, j = 1, 2, …, s, and s is the total number of extracted facades.

[0070] The step 4 is specifically:

[0071] Step 4.1: Facade projection fitting

[0072] In the spatial rectangular coordinate system, the general plane equation is written as Ax + By + Cz + D = 0, and the plane parallel to the XOZ coordinate plane is taken as the projection plane; the point cloud data P = {p0, p1, …, p n} of the building facade is projected, first the minimum value y min and the maximum value y max of the Y-axis of the point cloud data are determined, and y = (y min +y max ) / 2 is taken as the projection plane; then the point cloud data is projected to the plane to obtain the projected point coordinates; finally, the least squares method is used to estimate the position of the main wall surface;

[0073] Step 4.2: Point distance constraint, points outside the given threshold range are deleted according to the distance from the point to the straight line;

[0074] Step 4.2.1, for each point in the projected point set P' = {p'0, p'1, …, p' n}, the distance from the point to the straight line is calculated;

[0075] Step 4.2.2, if the distance is less than the distance threshold d w , the index of the point is stored in the array; otherwise, skip and continue to judge the next point until all points are judged;

[0076] Step 4.2.3, the main wall surface points are extracted from the original point cloud data P according to the point index, and the main wall surface detection of the building facade is realized.

[0077] Step 5 is specifically:

[0078] Step 5.1, preliminary positioning of the window area through spatial grid division

[0079] Step 5.1.1, according to the size of the bounding box of the main wall surface obtained in step 4, a 3D space that can accommodate all point cloud data is constructed, and the 3D space is evenly divided into M*N*L small spatial grids along the X, Y and Z axes, and each spatial grid is numbered;

[0080] Step 5.1.2, the grid length of the point cloud data along three directions is calculated according to the number of grid division respectively, and the point cloud data is filled into the space grid to obtain the grid index corresponding to each point. The grid index (i, j, k) corresponding to the point (x, y, z) is calculated by formula (6). After the calculation is completed, the points with the same grid index are placed in the same space grid;

[0081]

[0082] In formula (6): cell x The size of the space grid in the X-axis direction; cell y The size of the space grid in the Y-axis direction; cell z The size of the space grid in the Z-axis direction;

[0083] Step 5.1.3, traverse the space grid. When there is no point in the grid, it is considered that the grid is located in the potential window area, and it is marked as 1. Otherwise, it is considered to be a wall, and it is marked as 0.

[0084] Step 5.1.4, count the number of grids along X and Y directions in the potential window area respectively, calculate the ratio and area, if the ratio is less than the set threshold T r and the area is greater than the threshold T s , it is considered that the region is an interference element, and the grid in the region is marked as 0; otherwise, it is considered that the region is a window region, and the grid index of the four corners is recorded.

[0085] Step 5.1.5, repeat step 5.1.4 until all potential window areas are judged, so as to obtain the approximate position of the candidate window in the building facade;

[0086] Step 5.1.6, according to the obtained preliminary candidate window position, taking the window boundary as the starting point, expanding the range of one grid to obtain the candidate window boundary grid. For the points in the candidate window boundary grid, store their indexes in the array V to avoid global point cloud participation in detection.

[0087] Step 5.2, extract the accurate boundary of the window and perform single body operation;

[0088] Step 5.2.1, first extract the point cloud subset R according to the index array V obtained in step 5.1. Select two points P1 and P2 in the subset R, draw a circle with a radius of a passing through the points P1 and P2, and judge whether the remaining points in R are in the circle one by one. If there is no point in the circle, it is considered that the points P1 and P2 are boundary points. The radius a is calculated by formula (7) and formula (8);

[0089]

[0090]

[0091] In formula (7) and (8): k-number of neighboring points; d-distance from neighboring points to the current point i

[0092] Step 5.2.2, repeat step 5.2.1 to traverse all points in the subset R, and the boundary points of the opening in the building facade can be obtained;

[0093] Step 5.2.3, single the boundary points obtained in step 5.2.2, according to the boundary point set s in ={u0,u1,…,u n}, calculate its centroid, and find the point closest to the centroid by KNN algorithm, and take the point as the initial seed point;

[0094] Step 5.2.4, calculate the nearest Euclidean distance of the remaining points to the initial seed point according to formula (9), and then calculate the probability of the remaining points being selected as the next cluster center by formula (10), and take the point with the maximum probability as the next cluster center;

[0095] D(c i )=min(dist(c i ,u1),dist(c i ,u2),dist(c i ,u3),…,dist(c i ,u n ))(9)

[0096]

[0097] Step 5.2.5, repeatedly repeat step 5.2.4 until k cluster centers are obtained or the current distance change is less than the threshold value;

[0098] Step 5.2.6, calculate the distance between the remaining points and the k cluster centers, and divide the points closer to the center into the same class;

[0099] Step 5.2.7, recalculate the cluster centers according to the updated classification set;

[0100] Step 5.2.8, repeat steps 5.2.6 and 5.2.7 until the iteration number is reached or the cluster centers no longer change.

[0101] Step 5.2.9, remove the interference factors that do not belong to the window elements, set the area threshold γ S , the length-width ratio γ L , and the height threshold from the ground γ H , delete the interference elements, and realize the single extraction of the window.​

[0102] Step 6 is specifically:

[0103] Step 6.1, extracting feature points;

[0104] Step 6.1.1, extracting feature points by using QuickHull convex hull algorithm, first initializing a quadrilateral composed of four coordinate extreme points to form an initial convex hull;

[0105] Step 6.1.2, deleting the points in the boundary point set obtained in step 5.2 that are located inside the initial convex hull, and then finding the point farthest from the edge of the initial convex hull from the remaining points, taking the point as a new convex hull vertex, repeating the process until all points are processed, thus obtaining the required feature points;

[0106] Step 6.2, sorting and adjusting feature points;

[0107] Step 6.2.1, after obtaining the feature points in step 6.1.2, first calculating the centroid point O of the feature points;

[0108] Step 6.2.2, calculating the angle between the connecting line of each point to point O and the positive direction of the X axis, and sorting from large to small according to the size of the angle;

[0109] Step 6.2.3, according to the angle between the vector between adjacent two feature points and the X axis and the Y axis, it is judged whether the point is located in the horizontal direction or the vertical direction, when located in the horizontal direction, the y value of the point in the improper position needs to be placed in the y value of the point with the smallest angle with the X axis in the current direction; when located in the vertical direction, the x value of the point in the improper position needs to be placed in the x value of the point with the smallest angle with the Y axis in the current direction;

[0110] Step 6.2.4, in order to prevent arc-shaped windows from being transformed into rectangular windows, it is still necessary to calculate the ratio μ of the actual perimeter of each window to the ideal perimeter according to formula (11), if μ < 0.90, the window is an arc-shaped window, and no coordinate adjustment operation is performed on the feature points in the arc-shaped part of the window;

[0111]

[0112] In formula (11):

[0113] p i —A point in the feature points;

[0114] n—The number of feature points;

[0115] max y , min y —The maximum y value and the minimum y value in the monomer window point cloud data, respectively;

[0116] max x min x max and min are the maximum and minimum x values in the monomer window point cloud data, respectively;

[0117] Step 6.3, equidistant interpolation points;

[0118] Step 6.3.1, according to the arrangement order of the points in step 6.2.2, the distance between two adjacent points is calculated, and new points are inserted on the vector between the adjacent points according to the unit vector between the adjacent points and the distance between each new point, so as to obtain a refined boundary extraction, and finally the boundary points are connected two by two to generate a window contour line.

[0119] The beneficial effects of the present application are:

[0120] The window detection method based on hole constraint and hierarchical positioning provided by the present application first extracts the building facade in a complex outdoor scene; then, with the help of the hole idea, the main wall surface extracted from the building facade is subjected to hierarchical judgment to obtain accurate window boundary information, and the window boundary points are individualized through clustering; then, for each individual window boundary point, feature points are extracted, and the boundary points are regularized by equidistant interpolation points; finally, a window frame model is constructed to obtain the window detection result. The present application solves the problem of inaccurate window detection caused by irregular distribution and non-uniform shape of windows in a complex facade, and can realize window detection of complex buildings with different facade structures, thereby improving the detection accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0121] Figure 1 is the original data input in step 1 of the embodiment of the present application;

[0122] Figure 2 is the scene rotation result obtained in step 1 of the embodiment of the window detection method based on hole constraint and hierarchical positioning of the present application;

[0123] Figure 3 is the result graph of ground filtering in step 2 of the embodiment of the window detection method based on hole constraint and hierarchical positioning of the present application;

[0124] Figure 4 is the result graph of extracting buildings in step 3.1 of the embodiment of the window detection method based on hole constraint and hierarchical positioning of the present application;

[0125] Figure 5 is the result graph of extracting building facades in step 3.2 of the embodiment of the window detection method based on hole constraint and hierarchical positioning of the present application;

[0126] Figure 6is the main wall surface extraction result graph of step 4 of an embodiment of the window detection method of the application based on hole constraint and hierarchical positioning;

[0127] Figure 7 is the schematic diagram of step 5.1 of the preliminary window area positioning of an embodiment of the window detection method of the application based on hole constraint and hierarchical positioning;

[0128] Figure 8 is the window singulation result graph of step 5.2 of an embodiment of the window detection method of the application based on hole constraint and hierarchical positioning;

[0129] Figure 9 is the window boundary optimization result graph of step 6 of an embodiment of the window detection method of the application based on hole constraint and hierarchical positioning;

[0130] Figure 10 is the final window contour result graph of step 6 of an embodiment of the window detection method of the application based on hole constraint and hierarchical positioning. DETAILED DESCRIPTION

[0131] The application will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0132] Embodiment 1

[0133] The embodiment provides a window detection method based on hole constraint and hierarchical positioning, which is specifically implemented according to the following steps:

[0134] Step 1: filtering and rotating the point cloud original data as shown in Figure 1 , as shown in Figure 2 , so as to rotate the scene to be parallel to the ground and the horizontal plane;

[0135] Step 1.1: performing a hybrid filtering operation on the original point cloud data

[0136] Step 1.1.1: dividing the original point cloud data into a plurality of cubes with a side length of l, the value of the side length l being 0.02m-0.04m, calculating the center of gravity of the cube and using the center of gravity point to represent all points in the cube, and using the point in the original point cloud data closest to the center of gravity point to replace the center of gravity point;

[0137] Step 1.1.2: continuing to remove noise points in a radius filtering manner: judging whether the number of neighboring points in the specified radius sphere is less than a specified threshold as the basis for whether the point is a discrete point, if the number of neighboring points in the specified radius sphere is less than the specified threshold, the point is determined to be a discrete noise point and is deleted; otherwise, it is retained;

[0138] Step 1.2: realizing scene rotation

[0139] Step 1.2.1: Based on the feature that the normal vector of the ground point is parallel to the y-axis, the angle between the normal vector and the y-axis in the point cloud data is less than The points are preliminarily extracted, and a plane is fitted to obtain the normal vector n1 of the plane;

[0140] Step 1.2.2, given the normal vector n2 of the target plane, since the target plane needs to be parallel to the horizontal plane XOZ, select n2 = (0, 0, 1);

[0141] Step 1.2.3, calculate the rotation angle and rotation axis: Calculate the angle θ between n1 and n2 according to formula (1), that is, the rotation angle; the rotation axis is to solve the vector perpendicular to the main plane normal vector n1 and the target plane normal vector n2.

[0142]

[0143] Step 1.2.4, solve the rotation matrix R according to the rotation axis and the angle θ, and transform the original point cloud scene by the rotation matrix R.

[0144] In step 1.2.1 is 10°.

[0145] Step 2: Based on step 1, remove the ground points in the scene, such as Figure 3 As shown, non-ground points are clustered and building points are extracted by analyzing the features of non-ground points;

[0146] Step 2.1: Invert the point cloud data obtained in step 1, generate the initial cloth grid according to the set grid size, and determine the number of grid particles. The starting position of the cloth is usually above the highest point of the inverted point cloud;

[0147] Step 2.2: Project the point cloud data and grid particles to the same horizontal plane, search for the nearest neighbor of each grid particle, and record its elevation value IHV before projection;

[0148] Step 2.3: For each grid particle, if it can move, calculate its displacement caused by the external force and compare the position of the particle after displacement with the IHV of its corresponding nearest neighbor. If the height of the particle is less than or equal to the elevation value, the particle is considered to be an immovable point and its height is set to IHV.

[0149] Step 2.4, calculate the displacement distance of grid particles caused by internal factors;

[0150] Step 2.5: Repeat steps 2.3 and 2.4 until the maximum number of iterations is reached or the moving distance of all particles is less than the preset value, then stop the iteration;

[0151] Step 2.6, calculate the distance between the grid particle and its nearest neighbor: if less than the distance threshold dis g (dis g Set to 0.05), it is considered as a ground point, otherwise as a non-ground point. Thus, the preliminary ground point set G can be obtained;

[0152] Step 2.7, for the set G, calculate its average height and the normal vector of each point;

[0153] Step 2.8, for each point in the set G, if the height of the point is higher than the average height, calculate the angle between its normal vector and the z-axis, when the angle is less than the set threshold a, set a to 20°, then skip; otherwise, it is classified into the non-ground point set and deleted from the ground point set;

[0154] Step 2.9, through K nearest neighbor (KNN) search, the number of non-ground points belonging to the non-ground point set around the ground seed point, if the number of non-ground points is greater than 0.85 times the total number of nearest neighbors, then the point is classified into the non-ground point and deleted from the ground point set; otherwise, the non-ground points around the seed point are classified into the ground point;

[0155] Repeat steps 2.8 and 2.9 until all seed points have been traversed.

[0156] Step 3: extract the required building facade;

[0157] Step 3.1: extract the building using the method of fusion dimension feature and roughness;

[0158] Step 3.2: further extract the building facade using the improved RANSAC method;

[0159] Step 4: extract the main wall surface in the extracted building facade;

[0160] Step 5: on the main wall surface, extract the window boundary point using the hierarchical method;

[0161] Step 6: perform window boundary optimization operation.

[0162] Embodiment 2

[0163] The embodiment provides a window detection method based on hole constraint and hierarchical positioning, which is based on the embodiment 1

[0164] Step 3.1 is specifically as follows:

[0165] Step 3.1.01, input the non-ground point set Q = {q1, q2, …, q m}, m is the total number of non-ground points, and a kd tree is created to find the nearest neighbor points N of the search point q through KNN q ; Next, calculate N q The covariance eigenvalues ​​λ1, λ2, λ3 of the inner point cloud (λ1≥λ2≥λ3); finally, the surface index λ is calculated by formula (2) P , divergence index λ S and rod index λ L , which is collectively referred to as dimensional features, that is, calculating each point q j The dimension feature f qj ;

[0166]

[0167] Step 3.1.02, label all points in set Q, set label to 0, and initialize the label of the cluster segmentation block segLab to 1;

[0168] Step 3.1.03, traverse all points, if the label of the current seed point is 0, find the nearest related point q through the neighbor point search NN , and store it in array V q The associated points are filtered by the dimension threshold, and the associated points with different characteristics from the seed points are removed from V q , thereby obtaining a new associated point array V n ;

[0169] Step 3.1.04, traverse array V n For all points in the dataset, if there is a point with a label other than 0, find the smallest label value, compare it with segLab, and assign the smallest value of the two to minSegLab; if there is no point with a label other than 0, directly assign segLab to minSegLab.

[0170] Step 3.1.05, for array V n Traverse again. If the label of the current point is greater than minSegLab, traverse all points again. If the label of the point is the same as the label of the current point, reset the label value of the point to the value of minSegLab.

[0171] In step 3.1.06, add 1 to the value of segLab and repeat steps 3.1.3-3.1.5 until the values ​​of all point labels are not equal to the set value, thus completing the coarse segmentation of the entire point cloud;

[0172] Step 3.1.07, each segmentation block F obtained in step 3.1.6 i(i e 1,2,...,N), the average height of each segment is calculated by formula (3) and (4) respectively and roughness R s ;

[0173]

[0174]

[0175] In formula (3) and (4): — the average height of the segment; Z j — the height of each point in the segment; n — the number of points in the segment; R s — roughness; N — the number of segments;

[0176] Step 3.1.08, each segment is regarded as a whole, a classification label is added, and the first segment is taken as a seed surface, and the seed surface is used as a starting point for searching;

[0177] Step 3.1.09, all adjacent segments of the surface are found by spatial connectivity analysis. When the similarity threshold of the seed surface and the adjacent segment is different, it is considered that the two belong to different segment categories and cannot be merged, and the next adjacent segment is selected; otherwise, the label of the adjacent segment is reset to the label of the seed surface, and the next adjacent segment is judged. The similarity threshold is set according to the roughness of the point cloud surface to judge the similarity of different segments;

[0178] Step 3.1.10, it is judged whether all adjacent segments have been traversed, if there are adjacent segments that have not been traversed, step 3.1.9 is repeated until the merging of all adjacent segments is completed; otherwise, the segment that has not been accessed is taken as a seed surface;

[0179] Step 3.1.11, it is judged whether all segments have been accessed, if there are still segments that have not been accessed, steps 3.1.9-3.1.10 are repeated; otherwise, the merging is ended;

[0180] Step 3.1.12, after the merging is completed, the point cloud data is classified according to the classification label;

[0181] Step 3.1.13, the number of points in each class cluster after merging is judged, a maximum threshold max and a minimum threshold min are set, max and min are set to 10000000 and 100 respectively, and the class clusters within the range are retained;

[0182] Step 3.1.14, combined with different ground features, the retained class clusters are screened by using elevation and geometric features, and the class clusters belonging to building points are extracted, such as Figure 4As shown, the precise extraction of buildings is realized. Different ground features are as follows: building features are regular in shape, high in height, and have large point cloud density in projection area; pole features such as electric poles and street lamps have certain height and are perpendicular to the ground, have small projection area and large point density; vegetation features have certain height and have uniform point density distribution in projection range.

[0183] Step 3.2 is specifically:

[0184] Step 3.2.1, in the building point cloud data obtained in step 3.1, three points are randomly selected to form a plane, and the parameters of the plane model are calculated;

[0185] Step 3.2.2, the distance of all points in the data set to the plane is calculated, if within the distance threshold d W (set the distance threshold to 0.15), the point is classified as an inlier point, otherwise as an outlier point;

[0186] Step 3.2.3, repeat steps 3.2.1-3.2.2 until the iteration threshold is reached, then extract the plane with the largest number of inlier points;

[0187] Step 3.2.4, count the total number of points in the plane, when less than the threshold S num , delete it, otherwise project the extracted plane to the YOZ plane, so that all points form a vertical line segment, sort the projection points according to the y value from small to large, and use formula (5) to calculate the distance between adjacent two points;

[0188]

[0189] In formula (5), P i —point in the facade W j , j = 1, 2, …, s; K—number of points in the facade W j .

[0190] Step 3.2.5, determine whether there is a vertical gap in the facade by determining whether the correlation threshold between points is greater than 3 times the median distance threshold Dis y , wherein each gap can divide the data into two different point clusters, if there is a gap, the facade is divided again to realize the optimization of the extraction result; otherwise, skip;

[0191] Step 3.2.6, project the extracted plane to the XOZ plane, so that all points form a horizontal line segment, sort the projection points according to the x value from small to large, and use formula (5) to calculate the distance between adjacent two points;

[0192] Step 3.2.7, determine whether there is a vertical gap in the facade by determining whether the correlation threshold between points is greater than 3 times the median distance threshold Disx To detect whether there are horizontal gaps in the facade, each gap can be divided into two different point clusters. If there is a gap, the facade is split again to optimize the extraction result; otherwise, it is skipped;

[0193] Step 3.2.8, repeat steps 3.2.1-3.2.7 until all potential facades are extracted from the building point cloud, and then according to the main wall features, the main wall features include height, perpendicular to the XOZ surface features, such as Figure 5 As shown, the required building facades are extracted from all potential facades and added to the set W, W = {W1, W2, ... W s}, where W j Denoted as the j-th facade, j = 1, 2, …, s, where s is the total number of extracted facades.

[0194] Step 4 is as follows:

[0195] Step 4.1: Elevation Projection Fitting

[0196] In the spatial rectangular coordinate system, the plane equation is written as Ax+By+Cz+D=0, and the plane parallel to the XOZ coordinate plane is used as the projection plane; for the building facade point cloud data P={p0,p1,…,p n}For projection, first determine the minimum Y-axis value y of the point cloud data min and the maximum value y max , and let y=(y min +y max ) / 2 is the projection plane; then the point cloud data is projected onto the plane to obtain the coordinates of the projected points; finally, the least squares method is used to estimate the position of the main wall;

[0197] Step 4.2: Point distance constraint, delete points outside the given threshold range based on the distance from the point to the line;

[0198] Step 4.2.1, for the projection point set P′={p′0,p′1,…,p′ n}For each point in the equation, calculate the distance from the point to the line;

[0199] Step 4.2.2, if it is less than the distance threshold d w (d w If the value is between 0.1 and 0.2, the index of the point is stored in the array; otherwise, it is skipped and the next point is judged until all points are judged.

[0200] Step 4.2.3, extract the main wall points from the original point cloud data P according to the point index, such as Figure 6 As shown, the main wall detection of the building facade is realized.

[0201] Example 3

[0202] The embodiment provides a window detection method based on hole constraint and hierarchical positioning. On the basis of the embodiment 2, step 5 is specifically as follows:

[0203] Step 5.1, preliminary positioning of the window region through spatial grid division

[0204] Step 5.1.1, according to the size of the main wall surface bounding box obtained in step 4, a 3D space capable of containing all point cloud data is constructed, the 3D space is evenly divided into M*N*L small spatial grids along the X, Y and Z axes, and each spatial grid is numbered, wherein, because the points in the main wall surface change little in the Z axis direction, L is defined as 1, that is, the space is divided into M*N*1 small spatial grids;

[0205] Step 5.1.2, the grid lengths of the point cloud data along the three directions are calculated according to the grid division numbers respectively, and the point cloud data is filled into the spatial grid to obtain the grid index corresponding to each point, and formula (6) is used to calculate the grid index (i, j, k) corresponding to the point (x, y, z), after the calculation, the points with the same grid index are placed in the same spatial grid;

[0206]

[0207] In formula (6), cell x is the size of the spatial grid in the X axis direction; cell y is the size of the spatial grid in the Y axis direction; and cell z is the size of the spatial grid in the Z axis direction.

[0208] Step 5.1.3, the spatial grid is traversed, when there is no point in the grid, it is considered that the grid is located in a potential window region, and is marked as 1, otherwise, it is considered to be a wall, and is marked as 0;

[0209] Step 5.1.4, the grid numbers of the potential window region along the X and Y directions are counted respectively, the ratio and the area are calculated, if the ratio is less than a set threshold T r and the area is greater than a threshold T s (T r is 0.08, and T s is 40), it is considered that the region is an interference element, the grid in the region is marked as 0; otherwise, it is considered that the region is a window region, and the grid indexes located at the four corners are recorded;

[0210] Step 5.1.5, step 5.1.4 is repeated until all potential window regions are judged, so that the approximate position of the candidate window in the building facade is obtained;

[0211] Step 5.1.6, ifFigure 7 As shown, according to the obtained preliminary candidate window position, a candidate window boundary grid is obtained by expanding the window boundary by one grid in all directions, and for the points in the candidate window boundary grid, their indices are stored in an array V to avoid global point cloud participation in detection;

[0212] Step 5.2, extracting the accurate boundary of the window and performing a single operation thereon;

[0213] Step 5.2.1, first extracting a point cloud subset R according to the index array V obtained in step 5.1, and then selecting two points P1 and P2 in the subset R, drawing a circle with a radius of a passing through the points P1 and P2, and judging whether the remaining points in R are within the circle one by one. If there is no point within the circle, it is considered that the points P1 and P2 are boundary points. The radius a is calculated by formula (7) and formula (8);

[0214]

[0215]

[0216] In formula (7) and (8), k is the number of neighboring points; d i is the distance from the neighboring point to the current point;

[0217] Step 5.2.2, repeating step 5.2.1 to traverse all points in the subset R, and then the boundary points of the opening in the building facade can be obtained;

[0218] Step 5.2.3, performing a single operation on the boundary points obtained in step 5.2.2, calculating the centroid according to the boundary point set s in ={u0, u1, …, u n}, finding the point closest to the centroid by KNN algorithm, and taking the point as the initial seed point;

[0219] Step 5.2.4, calculating the nearest Euclidean distance of the remaining points to the initial seed point according to formula (9), and then calculating the probability of the remaining points being selected as the next clustering center by formula (10), and taking the point with the maximum probability as the next clustering center;

[0220] D(c i )=min(dist(c i ,u1),dist(c i ,u2),dist(c i ,u3),…,dist(c i ,u n )) (9)

[0221]

[0222] Step 5.2.5, repeat step 5.2.4 until k cluster centers are obtained or the current distance variation is less than a threshold value;

[0223] Step 5.2.6, calculate the distance between the remaining points and the k cluster centers, and divide the points closer to the center into the same class;

[0224] Step 5.2.7, recalculate the cluster centers according to the updated classification set;

[0225] Step 5.2.8, repeat steps 5.2.6 and 5.2.7 until the number of iterations is reached or the cluster centers no longer change.

[0226] Step 5.2.9, remove interference factors that do not belong to window elements, set an area threshold γ in the code S less than 0.9 or greater than 4, aspect ratio γ L greater than 3 (3:1), and set a height threshold γ from the ground H less than 0.5, delete the interference elements, as shown in Figure 8 to achieve the single extraction of the window.

[0227] Step 6 is specifically:

[0228] Step 6.1, extract feature points;

[0229] Step 6.1.1, use the QuickHull convex hull algorithm to extract feature points, first initialize a quadrilateral composed of four coordinate extreme points to form an initial convex hull;

[0230] Step 6.1.2, delete the points in the boundary point set obtained in step 5.2 that are inside the initial convex hull, and then find the point farthest from the edge of the initial convex hull from the remaining points, take this point as a new convex hull vertex, repeat the process until all points are processed, thus obtaining the required feature points;

[0231] Step 6.2, feature point sorting and adjustment;

[0232] Step 6.2.1, after obtaining the feature points of step 6.1.2, first calculate the centroid point O of the feature points;

[0233] Step 6.2.2, calculate the angle between the line connecting each point to point O and the positive direction of the X axis, and sort from large to small according to the angle;

[0234] Step 6.2.3, according to the angle between the vector of the two adjacent feature points and the X axis and the Y axis, it is determined whether the point is located in the horizontal direction or the vertical direction, when located in the horizontal direction, the y value of the point in the improper position needs to be placed in the y value of the point with the smallest angle with the X axis in the current direction; when located in the vertical direction, the x value of the point in the improper position needs to be placed in the x value of the point with the smallest angle with the Y axis in the current direction;

[0235] Step 6.2.4, in order to prevent the arc-shaped window from being transformed into a rectangular window, it is still necessary to calculate the ratio μ of the actual perimeter of each window to the ideal perimeter according to formula (11), if μ < 0.90, the window is an arc-shaped window, and the coordinate adjustment operation is not performed on the feature points in the arc-shaped part of the window;

[0236]

[0237] In formula (11), p i —A point in the feature points; n—the number of feature points; max y , min y —The maximum y value and the minimum y value in the monomer window point cloud data respectively; max x , min x —The maximum x value and the minimum x value in the monomer window point cloud data respectively;

[0238] Step 6.3, equidistant interpolation;

[0239] Step 6.3.1, according to the arrangement order of the points in step 6.2.2, the distance between two adjacent points is calculated, according to the unit of the vector between adjacent points and the distance (the size of the distance is set to 0.02) between each new point, a new point is inserted on the vector, so as to achieve the purpose of obtaining refined boundary extraction, and the boundary optimization result is as shown in Figure 9 Finally, the boundary points are connected two by two, and the window contour line is generated as shown in Figure 10 .

Claims

1. A window detection method based on hole constraints and hierarchical positioning, characterized in that: Please follow the steps below to implement it: Step 1: Filter and rotate the original point cloud data so that the scene is rotated to be parallel to the ground and the horizontal plane; Step 2: Based on step 1, remove the ground points in the scene, cluster the non-ground points, and extract the building points by analyzing the features of the non-ground points; Step 3: Extract the required building facade; Step 4: Extract the main wall from the extracted building facade; Step 5: On the main wall, use the hierarchical method to extract the window boundary points; Step 6: Perform window boundary optimization operation; The step 6 is specifically as follows: Step 6.1, extract feature points; Step 6.1.1: Use the QuickHull convex hull algorithm to extract feature points. First, initialize a quadrilateral consisting of four coordinate extreme points to form the initial convex hull. Step 6.1.2: Delete the points inside the initial convex hull from the boundary points obtained in step 5. Then, find the point farthest from the edge of the initial convex hull from the remaining points. Use this point as the new convex hull vertex. Repeat this process until all points have been processed. In this way, the required feature points can be obtained. Step 6.2, feature point sorting and adjustment; Step 6.2.1: After obtaining the feature points in step 6.1.2, first calculate the centroid O of the feature points; Step 6.2.2, calculate the angle between the line connecting each point to point O and the positive direction of the X-axis, and sort the points from largest to smallest according to the size of the angle; Step 6.2.3: Determine whether the point is in the horizontal or vertical direction based on the angle between the vector between two adjacent feature points and the X and Y axes. If the point is in the horizontal direction, the y value of the inappropriately positioned point needs to be set to the y value of the point in the current direction with the smallest angle with the X axis. If the point is in the vertical direction, the x value of the inappropriately positioned point needs to be set to the x value of the point in the current direction with the smallest angle with the Y axis. In step 6.2.4, to prevent curved windows from being transformed into rectangular windows, the ratio μ of the actual perimeter to the ideal perimeter of each window still needs to be calculated according to formula (11). If μ is less than 0.90, the window is a curved window, and the coordinate adjustment operation is not performed on the feature points of the curved part of the window. In formula (11): p i ——one of the feature points; n——the number of feature points; max y 、min y ——respectively the maximum y value and minimum y value in the single window point cloud data; max x 、min x ——respectively the maximum x value and minimum x value in the single window point cloud data; Step 6.3, equidistant interpolation points; In step 6.3.1, the distance between two adjacent points is calculated based on the order of the points in step 6.2.

2. New points are inserted into the vector based on the vector unit between the adjacent points and the spacing between each new point to achieve the purpose of obtaining refined boundary extraction. Finally, the boundary points are connected in pairs to generate the window outline.

2. The window detection method based on hole constraint and hierarchical positioning according to claim 1, characterized in that: Step 1 is as follows: Step 1.1: Perform hybrid filtering on the original point cloud data Step 1.1.1: Divide the original point cloud data into several cubes with a side length of l. Calculate the center of gravity of the cube and use it to represent all points in the cube. Replace the center of gravity with the point closest to the center of gravity in the original point cloud data. Step 1.1.2: Use radius filtering to continue removing noise points: determine whether the number of neighboring points in a sphere with a specified radius is less than a specified threshold as the basis for determining whether the point is a discrete point. If the number of neighboring points in a sphere with a specified radius is less than the specified threshold, the point is determined to be a discrete noise point and is deleted; otherwise, it is retained. Step 1.2: Implement scene rotation Step 1.2.1: Based on the feature that the normal vector of the ground point is parallel to the y-axis, the angle between the normal vector and the y-axis in the point cloud data is less than The points are preliminarily extracted, and a plane is fitted to obtain the normal vector n1 of the plane; Step 1.2.2, give the normal vector n2 of the target plane; Step 1.2.3, calculate the rotation angle and rotation axis: Calculate the angle θ between n1 and n2 according to formula (1), that is, the rotation angle; the rotation axis is to solve the vector perpendicular to the main plane normal vector n1 and the target plane normal vector n2. θ=arccos(n1,n2)(1); Step 1.2.4, solve the rotation matrix R according to the rotation axis and the angle θ, and transform the original point cloud scene by the rotation matrix R.

3. The window detection method based on hole constraints and hierarchical positioning according to claim 2, characterized in that: In step 1.2.1 is 10°.

4. The window detection method based on hole constraint and hierarchical positioning according to claim 1, characterized in that: The step 2 is specifically as follows: Step 2.1: Invert the point cloud data obtained in step 1, generate the initial cloth grid according to the set grid size, and determine the number of grid particles. The starting position of the cloth is usually above the highest point of the inverted point cloud; Step 2.2: Project the point cloud data and grid particles to the same horizontal plane, search for the nearest neighbor of each grid particle, and record its elevation value IHV before projection; Step 2.3: For each grid particle, if it can move, calculate its displacement caused by the external force and compare the position of the particle after displacement with the IHV of its corresponding nearest neighbor. If the height of the particle is less than or equal to the elevation value, the particle is considered to be an immovable point and its height is set to IHV. Step 2.4, calculate the displacement distance of grid particles caused by internal factors; Step 2.5: Repeat steps 2.3 and 2.4 until the maximum number of iterations is reached or the moving distance of all particles is less than the preset value, then stop the iteration; Step 2.6, calculate the distance between the grid particle and its nearest neighbor: if it is less than the distance threshold dis g , it is considered to be a ground point, otherwise it is a feature point, thus obtaining the preliminary ground point set G; Step 2.7, for the set G, calculate its average height and the normal vector of each point; Step 2.8: For each point in the set G, if the height of the point is higher than the average height, calculate the angle between its normal vector and the z-axis. If the angle is less than the set threshold α, skip it; otherwise, classify it as a non-ground point set and delete it from the ground point set. Step 2.9, search the number of non-ground points around the ground seed point through K-nearest neighbor (KNN); Repeat steps 2.8 and 2.9 until all seed points have been traversed.

5. The window detection method based on hole constraint and hierarchical positioning according to claim 1, characterized in that: The step 3 is specifically as follows: Step 3.1: Extract buildings by fusing dimensionality features and roughness. Step 3.2: Use the improved RANSAC method to further extract the building facade.

6. The window detection method based on hole constraint and hierarchical positioning according to claim 5, characterized in that: The step 3.1 is as follows: Step 3.1.01, input the non-ground point set Q obtained in step 2 = {q1,q2,…,q m }, m is the total number of non-ground points, and a kd tree is created to find the nearest neighbor points N of the search point q through KNN q ; Next, calculate N q The covariance eigenvalues ​​λ1, λ2, λ3 of the inner point cloud (λ1≥λ2≥λ3); finally, the surface index λ is calculated by formula (2) P , divergence index λ S and rod index λ L , which is collectively referred to as dimensional features, that is, calculating each point q j The dimension feature f qj ; Step 3.1.02, label all points in set Q, set label, and initialize the label of cluster segmentation block to segLab; Step 3.1.03, traverse all points, if the label of the current seed point is equal to the set value, find the nearest related point q through the neighbor point search NN , and store it in array V q In the process, the associated points are screened by the dimension threshold, and the associated points with different characteristics from the seed points are removed from V q , thereby obtaining a new associated point array V n ; Step 3.1.04, traverse array V n For all points in the dataset, if there is a point whose label is not equal to the set value, find the minimum label value, compare it with segLab, and assign the minimum value of the two to minSegLab; if there is no point whose label is not equal to the set value, directly assign segLab to minSegLab; Step 3.1.05, for array V n Traverse again. If the label of the current point is greater than minSegLab, traverse all points again. If the label of a point is the same as the label of the current point, reset the label value of the point to the value of minSegLab. In step 3.1.06, add 1 to the value of segLab and repeat steps 3.1.3-3.1.5 until the values ​​of all point labels are not equal to the set value, thus completing the coarse segmentation of the entire point cloud; Step 3.1.07, each segmentation block F obtained in step 3.1.6 i (i∈1,2,…,N), calculate the average height of each segmented block using formulas (3) and (4) respectively and roughness R s ; In formulas (3) and (4): ——the average height of the segmented blocks; Z j ——the height of each point in the segmentation block; n——the number of points in the segmentation block; R s ——Roughness; N——Number of segments; In step 3.1.08, each segment is considered as a whole, a classification label is added to it, and the first segment is used as the seed surface from which the query begins. Step 3.1.09: Use spatial connectivity analysis to find all adjacent segments of the face. If the similarity thresholds between the seed face and the neighboring segments are different, it can be assumed that the two segments belong to different segment categories and cannot be merged. In this case, the next adjacent segment is selected. Otherwise, the label of the adjacent segment is reset to the seed face label, and the next adjacent segment is judged. The similarity threshold is set according to the surface roughness of the point cloud to judge the similarity of different segments. Step 3.1.10: Determine whether all neighboring blocks have been traversed. If there are still neighboring blocks that have not been traversed, repeat step 3.1.9 until all neighboring blocks are merged. Otherwise, use the unvisited blocks as seed faces. Step 3.1.11: Check whether all the partitions have been visited. If there are still partitions that have not been visited, repeat steps 3.1.9-3.1.10; otherwise, the merge is terminated. Step 3.1.12: After the merging is completed, classify the point cloud data according to the classification labels; Step 3.1.13: Determine the number of points in each cluster after merging, set a maximum number threshold (max) and a minimum number threshold (min), and retain clusters within this range; In step 3.1.14, the retained clusters are screened using elevation and geometric features in combination with different landform features, and clusters belonging to building points are extracted to achieve accurate extraction of buildings.

7. The window detection method based on hole constraint and hierarchical positioning according to claim 5, characterized in that: Step 3.2 is as follows: Step 3.2.1, randomly select multiple points from the building point cloud data obtained in step 3.1 to form a plane, and calculate the parameters of the plane model; Step 3.2.2, calculate the distance from all points in the data set to the plane. If the distance threshold d W If the point is within the range, it is classified as an inside point, otherwise it is classified as an outside point. Step 3.2.3: Repeat steps 3.2.1-3.2.2 until the number of iterations reaches the iteration threshold, and then extract the plane with the largest number of local points; Step 3.2.4, count the total number of plane points, when it is less than the threshold S num If , it will be deleted. Otherwise, the extracted plane will be projected onto the YOZ plane so that all points form a vertical line segment. The projected points will be sorted from small to large according to the y value, and the distance between two adjacent points will be calculated using formula (5); In formula (5): P i ——Facade W j Points in, j = 1, 2, ..., s; K——Facade W j the number of midpoints; Step 3.2.5, by determining which points have a correlation threshold greater than 3 times the median distance threshold Dis y To detect whether there are vertical gaps in the facade, each gap can be divided into two different point clusters. If there is a gap, the facade is split twice to optimize the extraction results; otherwise, it is skipped. Step 3.2.6, project the extracted plane onto the XOZ plane so that all points form a horizontal line segment, sort the projected points in ascending order of x value, and use formula (5) to calculate the distance between two adjacent points; Step 3.2.7, by determining which points have a correlation threshold greater than 3 times the median distance threshold Dis x To detect whether there are horizontal gaps in the facade, each gap can be divided into two different point clusters. If there is a gap, the facade is split again to optimize the extraction results; otherwise, it is skipped; Step 3.2.8, repeat steps 3.2.1-3.2.7 until all potential facades are extracted from the building point cloud, and then extract the required building facade from all potential facades based on the main wall features, including height and perpendicularity to the XOZ plane features, and add them to the set W, W = {W1, W2, ... W s }, where W j Denoted as the j-th facade, j = 1, 2, …, s, where s is the total number of extracted facades.

8. The window detection method based on hole constraints and hierarchical positioning according to claim 1, characterized in that: The step 4 is specifically as follows: Step 4.1: Elevation Projection Fitting In the spatial rectangular coordinate system, the general form of the plane equation is written as Ax+By+Cz+D=0, and the plane parallel to the XOZ coordinate plane is used as the projection plane; For building facade point cloud data P = {p0,p1,…,p n }For projection, first determine the minimum Y-axis value y of the point cloud data min and the maximum value y max , and let y=(y min +y max ) / 2 is the projection plane; then the point cloud data is projected onto the plane to obtain the coordinates of the projected points; finally, the least squares method is used to estimate the position of the main wall; Step 4.2: Point distance constraint, delete points outside the given threshold range based on the distance from the point to the line; Step 4.2.1, for the projection point set P′={p0′,p1′,…,p n '}, calculate the distance from the point to the line; Step 4.2.2, if it is less than the distance threshold d w , then store the index of the point in the array; otherwise skip it and continue to judge the next point until all points are judged; In step 4.2.3, the main wall points are extracted from the original point cloud data P according to the point index to realize the main wall detection of the building facade.

9. The window detection method based on hole constraint and hierarchical positioning according to claim 1, characterized in that: The step 5 is specifically as follows: Step 5.1: Preliminary positioning of window areas through spatial grid division Step 5.1.1: Based on the main wall bounding box size obtained in step 4, construct a 3D space that can accommodate all point cloud data. Divide the 3D space into M*N*L small space grids along the X, Y, and Z axes, and number each space grid. Step 5.1.2: Calculate the grid side lengths of the point cloud data along the three directions according to the number of grid divisions, and fill the point cloud data into the spatial grid to obtain the grid index corresponding to each point. Use formula (6) to calculate the grid index (i, j, k) corresponding to the point (x, y, z). After the calculation is completed, points with the same grid index are placed in the same spatial grid. In formula (6): cell x ——The size of the spatial grid in the X-axis direction; cell y ——The size of the spatial grid in the Y-axis direction; cell z ——The size of the spatial grid in the Z-axis direction; Step 5.1.3, traverse the spatial grid. If there is no point in the grid, the grid is considered to be in the potential window area and marked as 1. Otherwise, it is considered to be a wall and marked as 0. Step 5.1.4: Count the number of grids in the potential window area along the X and Y directions, calculate their ratio and area, and if the ratio is less than the set threshold T r And the area is greater than the threshold T s , then the area is considered to be an interference element and the grid in the area is marked as 0; otherwise, the area is considered to be a window area and the grid indexes at the four corners are recorded; Step 5.1.5: Repeat step 5.1.4 until all potential window areas are determined, thereby obtaining the approximate locations of candidate windows on the building facade; Step 5.1.6: Based on the obtained preliminary candidate window positions, starting from the window boundary, expand the range of one grid in all directions to obtain the candidate window boundary grid. For each point in the candidate window boundary grid, store its index into array V to avoid global point cloud participation in detection. Step 5.2: extract the precise boundary of the window and perform a monomer operation on it; Step 5.2.1: First, extract the point cloud subset R according to the index array V obtained in step 5.

1. In the subset R, randomly select two points P1 and P2, draw a circle with a radius of α through points P1 and P2, and determine whether the remaining points in R are within the circle. If there is no point within the circle, then points P1 and P2 are considered to be boundary points, where the radius α is calculated by formula (7) and formula (8); In formulas (7) and (8): k is the number of neighboring points; d i ——The distance from the neighboring point to the current point; Step 5.2.2: Repeat step 5.2.1 to traverse all points in subset R to obtain the boundary points of the openings in the building facade. Step 5.2.3: Perform a monomer operation on the boundary points obtained in step 5.2.2, and calculate the boundary points according to the boundary point set s. in ={u0,u1,…,u n }, calculate its centroid, and use the KNN algorithm to find the point closest to the centroid and use this point as the initial seed point; Step 5.2.4, calculate the nearest Euclidean distance from the remaining points to the initial seed point according to formula (9), and then use formula (10) to calculate the probability of the remaining points being selected as the next cluster center, and select the point with the highest probability as the next cluster center; D(c i )=min(dist(c i ,u1),dist(c i ,u2),dist(c i ,u3),…,dist(c i ,u n ))(9) Step 5.2.5: Repeat step 5.2.4 until k cluster centers are obtained or the current distance change is less than the threshold; Step 5.2.6, calculate the distance between the remaining points and the k cluster centers, and classify the points closer to the center into the same class; Step 5.2.7, recalculate the cluster centers based on the updated classification set; Step 5.2.8: Repeat steps 5.2.6 and 5.2.7 until the number of iterations is reached or the cluster center no longer changes; Step 5.2.9, remove interference factors that are not window elements, and set the area threshold γ in the code S , aspect ratio γ L , and set the height threshold γ from the ground H , delete the interfering elements and realize the individual extraction of windows.

Images