Building surface parallax unsupervised optimization method based on local gradient descent

Abstract

The invention discloses a building surface parallax unsupervised optimization method based on local gradient descent, and belongs to the technical field of computer vision and urban satellite image three-dimensional reconstruction. The problem that parallax errors exist on the edges and the surfaces of buildings in urban satellite images is solved. According to the method, an urban satellite image and a parallax image are jointly used, color, spatial position and parallax information are coupled, and the color, spatial position and parallax information are mapped to a 6-dimensional feature space. And performing K-means clustering on building superpixels with obvious differences by using six-dimensional features. The superpixel areas are used as units, a local gradient descent algorithm is used for iteratively optimizing the parallax value in each superpixel area, the homography model corresponding to the superpixel areas is used as a target function in the optimization process, the parallax optimization process fully follows the geometric structure of a building, the distance between the parallax value of bilateral filtering and the homography model is used as a loss function, and the parallax value in each superpixel area is obtained. And reverse interference on optimization caused by an original parallax error is avoided. The method can be applied to building surface parallax optimization.

Description

Technical Field

[0001] The invention belongs to the technical field of computer vision and three-dimensional reconstruction of urban satellite images, and in particular relates to an unsupervised optimization method for building surface parallax based on local gradient descent. Background Art

[0002] In the field of aerospace photogrammetry research, 3D photogrammetry has become a research hotspot and future research trend, thanks to its advantages in wide coverage, high timeliness, high precision and multi-dimensional data acquisition. Therefore, aerospace 3D photogrammetry has become an important technical means for applications such as environmental monitoring and protection, disaster monitoring and emergency response, and urban planning and management. In addition, as the economic, cultural, political, technological and social service center of human society, cities play a vital role in the development and progress of society. The 3D reconstruction of cities based on aerospace photogrammetry can provide strong technical support and data guarantee for many fields such as urban planning, smart city construction, cultural heritage protection, education and scientific research.

[0003] The 3D reconstruction method based on urban satellite images generally includes the steps of satellite image stereo correction, stereo matching, 3D point cloud computing, triangular mesh generation and texture mapping. Among them, stereo matching can provide a disparity map that closely corresponds to the 3D point cloud data. Therefore, stereo matching provides the most direct source of 3D information for the entire 3D reconstruction process and directly affects the accuracy of the 3D model. Research on the 3D reconstruction method of urban satellite images is mostly focused on the stereo matching stage of urban satellite images.

[0004] The proposed urban satellite image stereo matching methods can be divided into two categories: one is the neural network-based stereo matching method, and the other is the traditional stereo matching method. Thanks to the rapid development of deep learning technology in the field of computer vision, the neural network-based stereo matching method has been well inspired. First, the stereo matching method based on the multi-scale convolutional neural network (Multi-Scale CNN) extracts multi-scale image features and calculates the matching cost. Later, the end-to-end stereo matching method was proposed. This method directly inputs the left and right view pairs and outputs the disparity map. The Geometry and Context Network (GC-Net) is a typical end-to-end convolutional neural network. GC-Net uses a three-dimensional convolutional network (3D CNN) for cost aggregation, directly outputs the disparity map by regression, and has achieved excellent results on multiple public datasets. In order to avoid the influence of texture inconsistency and illumination changes, the Guided Aggregation Network (GA-Net) was proposed. GA-Net enhances the cost aggregation process by combining the guided aggregation layer and the deep convolutional neural network, thereby generating a more accurate disparity map. In order to reduce the complexity of the network, the Hierarchical Neural Architecture Search (HNAS) method was proposed. The HNAS method uses a hierarchical neural architecture search strategy to automatically find the optimal network structure suitable for stereo matching tasks. Different from the traditional manually designed network, HNAS finds an efficient and high-performance network architecture through the definition and optimization of the search space, which significantly improves the accuracy and speed of stereo matching.

[0005] Although the neural network-based stereo matching method can achieve excellent results in many public standard data sets, it is difficult to apply to actual satellite images that lack standard disparity results and have diverse parameters due to its dependence on real standard data and the weak generalization ability of neural networks. In the actual task of satellite image 3D reconstruction, traditional stereo matching methods are more often used. Traditional stereo matching methods generally go through four processes: matching cost calculation, matching cost aggregation, disparity calculation, and disparity optimization. Since the local matching cost calculation method combined with the matching cost aggregation method based on dynamic programming can extract local features for matching cost calculation and has low time complexity. Therefore, this type of method can adapt to large-scale satellite images containing complex and changeable ground scenes. The most popular methods include Semi-Global Matching (SGM) method, More Global Matching (MGM) method, and Double Propagation Matching (DPM) method. These methods first calculate the local features of the image by fixed window or adaptive window, and calculate the matching cost between stereo images. Then, the matching costs within a certain range or on a specific path are aggregated through a dynamic programming algorithm, with the aim of optimizing the matching costs in the depth discontinuity areas at the edge of the building and the weak texture areas on the surface of the building.

[0006] Although the existing methods can obtain satisfactory results to a certain extent, this optimization algorithm based on dynamic programming algorithm is a globally consistent optimization method, lacks clear guidance of edges and regions, and is easy to introduce "star" and "stripe" errors. Therefore, in order to further optimize the parallax errors of building edges and surfaces in urban satellite images and get rid of the dependence on standard parallax results, it is necessary to propose a new building surface parallax optimization method. Summary of the invention

[0007] The purpose of the present invention is to solve the problem of parallax errors in building edges and surfaces in urban satellite images, and to propose an unsupervised optimization method for building surface parallax based on local gradient descent.

[0008] The technical solution adopted by the present invention to solve the above technical problems is:

[0009] An unsupervised optimization method for building surface parallax based on local gradient descent, the method specifically comprising the following steps:

[0010] Step 1: collect the left view and right view of the urban building from different perspectives, obtain a disparity map based on the left view and the right view, and perform superpixel segmentation on the pixel points in the disparity map based on the information in the satellite image of the urban building and the disparity map to obtain a superpixel segmentation result;

[0011] Step 2: Initialize the number of iterations q = 1;

[0012] Step 3, initialize the number of superpixels k=1;

[0013] Step 4: Calculate the parameters of the homography model H according to the disparity map corresponding to the kth superpixel;

[0014] Step 5: performing bilateral filtering on the disparity map corresponding to the k-th superpixel to obtain a disparity map after bilateral filtering;

[0015] Step 6: According to the disparity map after bilateral filtering, the matching point (x′, y′) corresponding to the pixel point (x, y) in the left view in the right view is re-determined, (x, y) and (x′, y′) are the re-determined matching point pairs, and the loss function corresponding to the pixel point (x, y) is calculated according to the homography model parameters and the re-determined matching point pairs;

[0016] Similarly, the loss function value corresponding to each pixel in the left view is obtained respectively;

[0017] Step 7: Determine whether k=K is satisfied;

[0018] If satisfied, proceed to step eight;

[0019] Otherwise, set k=k+1 and return to step 4;

[0020] Step 8: Determine whether the loss function values ​​corresponding to each pixel of the disparity map are less than the disparity threshold ε disp ;

[0021] If the loss function values ​​corresponding to each pixel in the disparity map are all less than the disparity threshold, the disparity map obtained in the previous iteration is the final disparity map;

[0022] If the loss function values ​​corresponding to the pixels in the disparity map are not all less than the disparity threshold, execute step nine;

[0023] Step nine, determine whether the number of iterations q is equal to the threshold Q;

[0024] If the number of iterations q is equal to the threshold Q, the x′ coordinate of the matching point (x′, y′) is updated according to the loss function value corresponding to each pixel point, and the final disparity map is obtained according to the updated coordinates;

[0025] If the number of iterations q is not equal to the threshold Q, the x′ coordinate of the matching point (x′, y′) is updated according to the loss function value corresponding to each pixel point, and the disparity map is updated according to the updated coordinates, and then q=q+1 is set, and the process returns to step 3.

[0026] The beneficial effects of the present invention are:

[0027] 1. The present invention uses urban satellite images and disparity images in combination, couples color information, spatial position information and disparity information, and maps them into a 6-dimensional feature space. K-means clustering using 6-dimensional features can accurately segment building superpixels with obvious differences in color, distance and structure. Moreover, the superpixel segmentation method of the present invention is also applicable to applications such as ground object classification, recognition and map vectorization based on satellite images.

[0028] 2. The present invention uses the local gradient descent algorithm to iteratively optimize the disparity value in each superpixel region, taking the superpixel region as the unit. The process uses the homography model corresponding to the superpixel region as the objective function. The disparity optimization process fully follows the geometric structure of the building, so it can effectively optimize the disparity error in the depth discontinuity area of ​​the building edge, the disparity error in the weak texture area on the building surface, and the disparity hole error; the distance between the disparity value of the bilateral filter and the homography model is used as the loss function, avoiding the reverse interference of the original disparity error on the optimization.

[0029] 3. The building surface parallax optimization method based on local gradient descent proposed in the present invention is an unsupervised learning method, so it does not rely on standard parallax results and manual annotation results, saving a lot of manpower costs. Moreover, the parallax optimization method proposed in the present invention does not rely on a specific stereo matching method, so it is suitable for post-processing of various stereo matching algorithms. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 It is a flow chart of an unsupervised optimization method of building surface parallax based on local gradient descent of the present invention;

[0031] Figure 2 It is a flow chart of superpixel segmentation of coupled satellite image and disparity map;

[0032] Figure 3(a) is a satellite image;

[0033] Figure 3(b) is the disparity map corresponding to the satellite image;

[0034] Figure 3(c) is the superpixel segmentation result on the disparity map;

[0035] Among them, the red lines indicate the edges of superpixels.

[0036] Figure 4It is a flowchart of iterative optimization of building surface parallax based on local gradient descent;

[0037] Figure 5 This is a schematic diagram of iterative update of building surface parallax based on local gradient descent;

[0038] Among them, the blue thick line represents the objective function, the red curve represents the disparity after bilateral filtering, the green arrow represents the loss function gradient, the arrow length represents the learning rate, and the red dot represents the disparity value;

[0039] Figure 6(a) is a satellite image of the first test area;

[0040] FIG6( b ) is a disparity map of the first test area obtained by using the semi-global stereo matching algorithm;

[0041] FIG6( c ) is a disparity map of the first test area obtained by using the method of the present invention;

[0042] Figure 7(a) is a satellite image of the second test area;

[0043] FIG7( b ) is a disparity map of the second test area obtained by using the semi-global stereo matching algorithm;

[0044] FIG7( c ) is a disparity map of the second test area obtained by using the method of the present invention;

[0045] Figure 8(a) is a satellite image of the third test area;

[0046] FIG8( b ) is a disparity map of the third test area obtained by using the semi-global stereo matching algorithm;

[0047] FIG8( c ) is a disparity map of the third test area obtained by the method of the present invention; DETAILED DESCRIPTION

[0048] Specific implementation method 1: Combination Figure 1 The present embodiment describes an unsupervised optimization method for building surface parallax based on local gradient descent, and the method specifically comprises the following steps:

[0049] Step 1: collect the left view and right view of the urban building from different perspectives, and then obtain a disparity map based on the left view and the right view (the disparity map is obtained by a conventional stereo matching method), and perform super-pixel segmentation on the pixels in the disparity map based on the information in the satellite image of the urban building and the disparity map to obtain a super-pixel segmentation result;

[0050] Step 2: Initialize the number of iterations q = 1;

[0051] Step 3, initialize the number of superpixels k=1;

[0052] Step 4: Calculate the parameters of the homography model H according to the disparity map corresponding to the kth superpixel;

[0053] Step 5: performing bilateral filtering on the disparity map corresponding to the k-th superpixel to obtain a disparity map after bilateral filtering;

[0054] Step 6: Re-determine the matching point (x′, y′) corresponding to the pixel point (x, y) in the left view in the right view according to the disparity map after bilateral filtering. Since the disparity only exists in the X direction, y=y′ in the matching point represented by the disparity. (x, y) and (x′, y′) are the re-determined matching point pairs. Calculate the loss function corresponding to the pixel point (x, y) according to the homography model parameters and the re-determined matching point pairs.

[0055] Similarly, the loss function value corresponding to each pixel in the left view is obtained respectively;

[0056] Step 7: Determine whether k=K is satisfied;

[0057] If satisfied, proceed to step eight;

[0058] Otherwise, set k=k+1 and return to step 4;

[0059] Step 8: Determine whether the loss function values ​​corresponding to each pixel of the disparity map are less than the disparity threshold ε disp , ε disp =2;

[0060] If the loss function values ​​corresponding to each pixel in the disparity map are all less than the disparity threshold, the disparity map obtained in the previous iteration is the final disparity map;

[0061] If the loss function values ​​corresponding to the pixels in the disparity map are not all less than the disparity threshold, execute step nine;

[0062] Step nine, determine whether the number of iterations q is equal to the threshold Q;

[0063] If the number of iterations q is equal to the threshold Q, the x′ coordinate of the matching point (x′, y′) is updated according to the loss function value corresponding to each pixel point, and the final disparity map is obtained according to the updated coordinates;

[0064] If the number of iterations q is not equal to the threshold Q, the x′ coordinate of the matching point (x′, y′) is updated according to the loss function value corresponding to each pixel point, and the disparity map is updated according to the updated coordinates, and then q=q+1 is set, and the process returns to step 3.

[0065] The present invention uses satellite images and original disparity maps to perform superpixel segmentation to divide buildings with different structures, and then establishes a local gradient descent optimization mechanism to iteratively optimize the disparity of different surface structures.

[0066] Specific implementation method 2: Combination Figure 2 The present embodiment is different from the first embodiment in that, in step 1, the pixel points in the disparity map are segmented into superpixels according to the information in the satellite image of the urban building and the disparity map to obtain the superpixel segmentation result; the specific process is as follows:

[0067] Step 1: Initialize the number of superpixels to K; Initialize the label matrix Met label The same size as the disparity map, each element value in the label matrix represents the superpixel to which the pixel at the corresponding position in the disparity map belongs, and the label matrix Met label The initial value of each element in is -1;

[0068] Initialize the distance matrix Met distance The same size as the disparity map, the distance matrix Met distance Each element value in represents the distance from the pixel point at the corresponding position in the disparity map to the seed point of the superpixel to which it belongs, and the distance matrix Met distance The initial value of each element in is infinite, which means that the initial distance between all pixels and the seed point is very large;

[0069] Then the initial size S of each superpixel is:

[0070]

[0071] Where N is the total number of pixels in the disparity map;

[0072] Step 12: Divide the entire disparity map into K super-pixel regions of size S, and the width and height of each super-pixel region are equal;

[0073] For example, for a disparity map with a width of W and a height of H, the pixels in the 1st column to the W / Kth column form the first superpixel region, the pixels in the W / K+1th column to the 2W / Kth column form the second superpixel region, and so on;

[0074] Step 13: Map each pixel in the disparity map to a 6-dimensional vector. For pixel i, map pixel i to a vector (x i ,y i ,L i ,a i ,b i ,d i );

[0075] Among them, x i ,y i represents the plane space coordinates of pixel i, L i ,ai ,b i represents the color coordinates of pixel i in the satellite image mapped to the CIELab color space, d i Represents the disparity value of pixel i in the disparity map;

[0076] Step 14: For any superpixel region, calculate the mean of the vectors corresponding to all pixels in the superpixel region to obtain the mean vector, and use the mean vector as the vector corresponding to the initial seed point of the superpixel region (x cen ,y cen ,L cen ,a cen ,b cen ,d cen );

[0077] in, S 2 is the total number of pixels in the superpixel area;

[0078] Similarly, the vectors corresponding to the initial seed points of each superpixel region are obtained respectively, that is, the initial seed points of each superpixel region are evenly distributed;

[0079] Step 15, initialize the number of iterations l = 1;

[0080] Step 16: for pixel point i in the disparity map, determine whether pixel point i exists in the 2S×2S area around each seed point, determine all seed points in the surrounding 2S×2S area where pixel point i exists, and then calculate the comprehensive distance between the vector corresponding to pixel point i and the vector corresponding to each determined seed point, and obtain the minimum comprehensive distance;

[0081] If the minimum comprehensive distance obtained is less than the pixel i in the distance matrix Met distance (x, y), then use the obtained minimum comprehensive distance to update the pixel i in the distance matrix Met distance The corresponding distance in (x, y) is calculated, and the pixel i is assigned to the superpixel corresponding to the minimum comprehensive distance. The label matrix Met is updated according to the superpixel to which the pixel i belongs. label The pixel value corresponding to the pixel i in ;

[0082] Similarly, each pixel in the disparity map is processed separately;

[0083] Step 17: For any superpixel, calculate the mean of the vectors corresponding to all the pixels belonging to the superpixel after the current iteration (i.e., calculate the mean in the plane space, color space, and disparity space respectively), and use the calculated mean vector as the vector corresponding to the new seed point of the superpixel;

[0084] Similarly, the vector corresponding to the new seed point of each superpixel is calculated respectively;

[0085] Step 18: Determine whether the number of iterations l is greater than 20, and for each superpixel, calculate the comprehensive distance between the vector corresponding to the new seed point obtained in the current iteration and the vector corresponding to the seed point obtained in the previous iteration;

[0086] If the number of iterations is greater than 20 or the comprehensive distances corresponding to the seed points of all superpixels are less than the threshold ε sup (ε sup =2), the final superpixel segmentation result is obtained, and the final label matrix Met is returned label ; As shown in FIG3(a) and FIG3(b), the satellite image and the disparity map are shown respectively, and the superpixel segmentation result obtained by the method of the present invention is shown in FIG3(c), where the label of each pixel indicates the superpixel region to which it belongs;

[0087] Otherwise, set the number of iterations l=l+1, and use the vectors corresponding to the new seed points of each superpixel to return to step 16.

[0088] The other steps and parameters are the same as those in the first embodiment.

[0089] The present invention couples the color information, plane space information and disparity information of the urban satellite image and maps them to a 6-dimensional feature space. Then, by means of K-means clustering, the seed point coordinates and pixel labels are iteratively updated until the seed point coordinates no longer change significantly or the specified number of iterations is reached, and the clustering results of different pixel points are obtained. The clustering results can divide the building areas with obvious differences in color, distance and structure, thereby obtaining the final superpixel segmentation results.

[0090] Specific implementation method three: This implementation method is different from specific implementation method one or two in that the comprehensive distance between the vector corresponding to the pixel point i and the vector corresponding to each determined seed point is calculated respectively, specifically:

[0091]

[0092] in, is the color space distance, is the plane space distance, is the parallax space distance;

[0093] Through and Perform weighted summation to calculate the comprehensive distance D i :

[0094]

[0095] Among them, α, β and γ represent and , and α+β+γ=1.

[0096] The other steps and parameters are the same as those in the first or second embodiment.

[0097] Specific implementation method four: Combination Figure 4 This embodiment is different from the first to third embodiments in that the specific process of step 4 is as follows:

[0098] Step 41: In the local gradient descent optimization method, the homography model H is used to characterize the building surface model. The homography model H is specifically:

[0099]

[0100] Among them, h i′j′ Represents the parameters in the homography model H, usually h 33 =1;

[0101] Step 42: Since the homography model can represent the mapping relationship between any planes, the matching point corresponding to the disparity map pixel i on the surface of a single building structure is will satisfy the mapping relationship of the same homography model, then the matching points (x i ,y i ) and (x i ′,y i ′) is expressed as:

[0102]

[0103] After finishing, we get:

[0104]

[0105] Rearrange the equation system of formula (8) to obtain the linear equation system:

[0106]

[0107] Then, the objective function H(θ) is established according to the linear equation system of formula (9):

[0108] H(θ)=argMin{err1,err2,...,err N} (10)

[0109] Where n = 1, 2, ..., N, N is the number of randomly selected matching point pairs, each group includes 4 matching point pairs, and the 4 matching point pairs are non-collinear (that is, at most 3 of the 4 matching point pairs are collinear), N = 4S2 , the intermediate variable err n The calculation based on the randomly selected nth group of matching points is as follows:

[0110]

[0111] The parameters of the homography model are obtained by solving the objective function.

[0112] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0113] In order to avoid the interference of error points, the present invention calculates the homography model parameters that best represent the building structure in the current superpixel area by minimizing the errors of multiple iterative fitting homography models.

[0114] Specific implementation mode 5: This implementation mode is different from the specific implementation modes 1 to 4 in that the specific process of the bilateral filtering is as follows:

[0115]

[0116] Among them, BF[d i ] represents the disparity value d of pixel i in the disparity map i The result of bilateral filtering, S represents the neighborhood centered on pixel i (usually a 7×7 neighborhood near the center point, but not limited to a 7×7 neighborhood), d i and d j Represent the disparity values ​​of pixel i and pixel j respectively, represents the weight of the spatial domain, represents the weight of the disparity domain, represents the normalized weight, ||·|| represents the 2-norm, and |·| represents the absolute value.

[0117] The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

[0118] In order to overcome the interference of disparity error on the objective function, the bilateral filtering function BF[·] is used to optimize the original disparity value. The bilateral filtering result combines the dual effects of the spatial domain and the disparity range domain. In flat areas, the disparity changes slightly, and the corresponding disparity range domain weight is close to 1, so the spatial domain weight plays a major role, approximating Gaussian filtering. At the edge of the building, the disparity changes greatly, and the disparity range domain weight becomes larger, thereby maintaining the edge information.

[0119] Specific implementation method 6: This implementation method is different from the specific implementation methods 1 to 5 in that the normalized weight for:

[0120]

[0121] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5.

[0122] Specific implementation method 7: This implementation method is different from any one of the specific implementation methods 1 to 6 in that the spatial domain weight The calculation method is as follows:

[0123]

[0124] Among them, (x, y) is the coordinate of pixel i, (u, v) is the coordinate of pixel j, σ s Represents the standard deviation of the spatial domain.

[0125] The other steps and parameters are the same as those in Specific Embodiments 1 to 6.

[0126] Specific implementation example 8: This implementation example is different from any one of specific implementation examples 1 to 7 in that the disparity domain weight The calculation method is as follows:

[0127]

[0128] Among them, σ r Represents the standard deviation of the disparity domain.

[0129] The other steps and parameters are the same as those in Specific Embodiments 1 to 7.

[0130] Specific implementation method 9: This implementation method is different from any one of specific implementation methods 1 to 8 in that the loss function corresponding to the pixel point (x, y) is calculated based on the parameters of the homography model and the re-determined matching point pair, specifically:

[0131]

[0132] Among them, LF(x,y) is the loss function corresponding to the pixel point (x,y).

[0133] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0134] Specific implementation method ten: This implementation method is different from any one of specific implementation methods one to nine in that the x' coordinate of the matching point (x', y') is updated according to the loss function value corresponding to each pixel point, and the disparity map is updated according to the updated coordinate; the specific process is:

[0135] Calculate the gradient of the loss function ▽LF(x,y):

[0136]

[0137] Then the x′ coordinate is updated as follows:

[0138] x″=x′-α·▽LF(x,y) (18)

[0139] Among them, x″ represents the update result of x′, α represents the learning rate, α=0.5, which is used to control the step size of each update. By updating the parameters in the opposite direction of the gradient to reduce the value of the loss function, the process diagram is as follows Figure 5 As shown;

[0140] Then in the updated disparity map, the disparity corresponding to the pixel point (x, y) is:

[0141] disp v ' al =disp val -α·▽LF(x,y) (19)

[0142] Among them, disp v ' al Indicates the disparity corresponding to the updated pixel point (x, y), disp val =x′-x.

[0143] The other steps and parameters are the same as those in Specific Embodiments 1 to 9.

[0144] This embodiment iteratively updates the disparity along the direction of the gradient of the loss function.

[0145] Example

[0146] This embodiment proposes an unsupervised optimization method for building surface parallax based on local gradient descent, and the method specifically includes the following steps:

[0147] Step 1: Collect the left view and right view of the urban building from different perspectives, and then obtain a disparity map based on the left view and the right view (the disparity map is obtained by a conventional stereo matching method);

[0148] And perform super-pixel segmentation on the pixel points in the disparity map according to the information in the satellite image of the urban building and the disparity map to obtain the super-pixel segmentation result;

[0149] The present invention couples the color information, plane space information and disparity information of the urban satellite image and maps them to a 6-dimensional feature space. Then, by means of K-means clustering, the seed point coordinates and pixel labels are iteratively updated until the seed point coordinates no longer change significantly or the specified number of iterations is reached, and the clustering results of different pixel points are obtained. The clustering results can divide the building areas with obvious differences in color, distance and structure, thereby obtaining the final superpixel segmentation results. The specific process is:

[0150] Step 1: Initialize the number of superpixels to K; Initialize the label matrix Metlabel The same size as the disparity map, each element value in the label matrix represents the superpixel to which the pixel at the corresponding position in the disparity map belongs, and the label matrix Met label The initial value of each element in is -1;

[0151] Initialize the distance matrix Met distance The same size as the disparity map, the distance matrix Met distance Each element value in represents the distance from the pixel point at the corresponding position in the disparity map to the seed point of the superpixel to which it belongs, and the distance matrix Met distance The initial value of each element in is infinite, which means that the initial distance between all pixels and the seed point is very large;

[0152] Then the initial size S of each superpixel is:

[0153]

[0154] Where N is the total number of pixels in the disparity map;

[0155] Step 12: Divide the entire disparity map into K super-pixel regions of size S, and the width and height of each super-pixel region are equal;

[0156] For example, for a disparity map with a width of W and a height of H, the pixels in the 1st column to the W / Kth column form the first superpixel region, the pixels in the W / K+1th column to the 2W / Kth column form the second superpixel region, and so on;

[0157] Step 13: Map each pixel in the disparity map to a 6-dimensional vector. For pixel i, map pixel i to a vector (x i ,y i ,L i ,a i ,b i ,d i );

[0158] Among them, x i ,y i represents the plane space coordinates of pixel i, L i ,a i ,b i represents the color coordinates of pixel i in the satellite image mapped to the CIELab color space, d i Represents the disparity value of pixel i in the disparity map;

[0159] Step 14: For any superpixel region, calculate the mean of the vectors corresponding to all pixels in the superpixel region to obtain the mean vector, and use the mean vector as the vector corresponding to the initial seed point of the superpixel region (x cen ,y cen ,L cen ,a cen ,b cen ,d cen );

[0160] in, S 2 is the total number of pixels in the superpixel area;

[0161] Similarly, the vectors corresponding to the initial seed points of each superpixel region are obtained respectively, that is, the initial seed points of each superpixel region are evenly distributed;

[0162] Step 15, initialize the number of iterations l = 1;

[0163] Step 16: for pixel point i in the disparity map, determine whether pixel point i exists in the 2S×2S area around each seed point, determine all seed points in the surrounding 2S×2S area where pixel point i exists, and then calculate the comprehensive distance between the vector corresponding to pixel point i and the vector corresponding to each determined seed point, and obtain the minimum comprehensive distance;

[0164] The comprehensive distance between the vector corresponding to the pixel point i and the vector corresponding to each determined seed point is calculated respectively, specifically:

[0165] Take any seed point as an example

[0166]

[0167] in, is the color space distance, is the plane space distance, is the parallax space distance;

[0168] Through and Perform weighted summation to calculate the comprehensive distance D i :

[0169]

[0170] Among them, α, β and γ represent and , and α+β+γ=1.

[0171] If the minimum comprehensive distance obtained is less than the pixel i in the distance matrix Met distance(x, y), then use the obtained minimum comprehensive distance to update the pixel i in the distance matrix Met distance The corresponding distance in (x, y) is calculated, and the pixel i is assigned to the superpixel corresponding to the minimum comprehensive distance. The label matrix Met is updated according to the superpixel to which the pixel i belongs. label The pixel value corresponding to the pixel i in ;

[0172] Similarly, each pixel in the disparity map is processed separately;

[0173] Step 17: For any superpixel, calculate the mean of the vectors corresponding to all the pixels belonging to the superpixel after the current iteration (i.e., calculate the mean in the plane space, color space, and disparity space respectively), and use the calculated mean vector as the vector corresponding to the new seed point of the superpixel;

[0174] Similarly, the vector corresponding to the new seed point of each superpixel is calculated respectively;

[0175] Step 18: Determine whether the number of iterations l is greater than 20, and for each superpixel, calculate the comprehensive distance between the vector corresponding to the new seed point obtained in the current iteration and the vector corresponding to the seed point obtained in the previous iteration;

[0176] If the number of iterations is greater than 20 or the comprehensive distances corresponding to the seed points of all superpixels are less than the threshold ε sup (ε sup =2), the final superpixel segmentation result is obtained, and the final label matrix Met is returned label ; As shown in FIG3(a) and FIG3(b), the satellite image and the disparity map are shown respectively, and the superpixel segmentation result obtained by the method of the present invention is shown in FIG3(c), where the label of each pixel indicates the superpixel region to which it belongs;

[0177] Otherwise, set the number of iterations l=l+1, and use the vectors corresponding to the new seed points of each superpixel to return to step 16.

[0178] The present invention uses satellite images and original disparity maps to perform superpixel segmentation to divide buildings with different structures, and then establishes a local gradient descent optimization mechanism to iteratively optimize the disparity of different surface structures.

[0179] Step 2: Initialize the number of iterations q = 1;

[0180] Step 3, initialize the number of superpixels k=1;

[0181] Step 4: Calculate the parameters of the homography model H according to the disparity map corresponding to the kth superpixel; specifically:

[0182] Step 41: In the local gradient descent optimization method, the homography model H is used to characterize the building surface model. The homography model H is specifically:

[0183]

[0184] Among them, h i′j′ Represents the parameters in the homography model H, usually h 33 =1;

[0185] Step 42: Since the homography model can represent the mapping relationship between any planes, the matching point corresponding to the disparity map pixel i on the surface of a single building structure is will satisfy the mapping relationship of the same homography model, then the matching points (x i ,y i ) and (x i ′,y i ′) is expressed as:

[0186]

[0187] After finishing, we get:

[0188]

[0189] Rearrange the equation system of formula (8) to obtain the linear equation system:

[0190]

[0191] Then, the objective function H(θ) is established according to the linear equation system of formula (9):

[0192] H(θ)=argMin{err1,err2,...,err N} (10)

[0193] Where n = 1, 2, ..., N, N is the number of randomly selected matching point pairs, each group includes 4 matching point pairs, and the 4 matching point pairs are non-collinear (that is, at most 3 of the 4 matching point pairs are collinear), N = 4S 2 , the intermediate variable err n The calculation based on the randomly selected nth group of matching points is as follows:

[0194]

[0195] The parameters of the homography model are obtained by solving the objective function.

[0196] Step 5: performing bilateral filtering on the disparity map corresponding to the k-th superpixel to obtain a disparity map after bilateral filtering;

[0197] In order to overcome the interference of disparity error on the objective function, the bilateral filtering function BF[·] is used to optimize the original disparity value. The bilateral filtering result combines the dual effects of the spatial domain and the disparity range domain. In flat areas, the disparity changes slightly, and the corresponding disparity range domain weight is close to 1, so the spatial domain weight plays a major role, approximating Gaussian filtering. At the edge of the building, the disparity changes greatly, and the disparity range domain weight becomes larger, thereby maintaining the edge information. The specific process of bilateral filtering is as follows:

[0198]

[0199] Among them, BF[d i ] represents the disparity value d of pixel i in the disparity map i The result of bilateral filtering, S represents the neighborhood centered on pixel i (usually a 7×7 neighborhood near the center point, but not limited to a 7×7 neighborhood), d i and d j Represent the disparity values ​​of pixel i and pixel j respectively, represents the weight of the spatial domain, represents the weight of the disparity domain, represents the normalized weight, ||·|| represents the 2-norm, and |·| represents the absolute value.

[0200]

[0201] Among them, (x, y) is the coordinate of pixel i, (u, v) is the coordinate of pixel j, σ s represents the standard deviation in the spatial domain, σ r Represents the standard deviation of the disparity domain.

[0202] Step 6: Based on the disparity map after bilateral filtering, re-determine the matching point (x′, y′) corresponding to the pixel point (x, y) in the left view in the right view. Since the disparity only exists in the X direction, y=y′ in the matching point represented by the disparity. (x, y) and (x′, y′) are the re-determined matching point pairs. Calculate the loss function corresponding to the pixel point (x, y) based on the homography model parameters and the re-determined matching point pairs:

[0203]

[0204] Among them, LF(x,y) is the loss function corresponding to the pixel point (x,y).

[0205] Similarly, the loss function value corresponding to each pixel in the left view is obtained respectively;

[0206] Step 7: Determine whether k=K is satisfied;

[0207] If satisfied, proceed to step eight;

[0208] Otherwise, set k=k+1 and return to step 4;

[0209] Step 8: Determine whether the loss function values ​​corresponding to each pixel of the disparity map are less than the disparity threshold ε disp , ε disp =2;

[0210] If the loss function values ​​corresponding to each pixel in the disparity map are all less than the disparity threshold, the disparity map obtained in the previous iteration is the final disparity map;

[0211] If the loss function values ​​corresponding to the pixels in the disparity map are not all less than the disparity threshold, execute step nine;

[0212] Step nine, determine whether the number of iterations q is equal to the threshold Q;

[0213] If the number of iterations q is equal to the threshold Q, the x′ coordinate of the matching point (x′, y′) is updated according to the loss function value corresponding to each pixel point, and the final disparity map is obtained according to the updated coordinates;

[0214] If the number of iterations q is not equal to the threshold Q, the x′ coordinate of the matching point (x′, y′) is updated according to the loss function value corresponding to each pixel point, and the disparity map is updated according to the updated coordinates, and then q=q+1 is set, and the process returns to step 3.

[0215] The specific process of updating the disparity map is:

[0216] Calculate the gradient of the loss function ▽LF(x,y):

[0217]

[0218] Then the x′ coordinate is updated as follows:

[0219] x″=x′-α·▽LF(x,y) (18)

[0220] Among them, x″ represents the update result of x′, α represents the learning rate, α=0.5, which is used to control the step size of each update. By updating the parameters in the opposite direction of the gradient to reduce the value of the loss function, the process diagram is as follows Figure 5 As shown;

[0221] Then in the updated disparity map, the disparity corresponding to the pixel point (x, y) is:

[0222] disp v ' al =disp val -α·▽LF(x,y) (19)

[0223] Among them, disp v ' al Indicates the disparity corresponding to the updated pixel point (x, y), disp val =x′-x.

[0224] Experimental Section

[0225] The satellite images of the three test areas are shown in Figures 6(a), 7(a) and 8(a), respectively. The disparity map corresponding to Figure 6(a) obtained by using the semi-global stereo matching algorithm is shown in Figure 6(b), and the disparity map corresponding to Figure 6(a) obtained by using the method of the present invention is shown in Figure 6(c). The disparity map corresponding to Figure 7(a) obtained by using the semi-global stereo matching algorithm is shown in Figure 7(b), and the disparity map corresponding to Figure 7(a) obtained by using the method of the present invention is shown in Figure 7(c). The disparity map corresponding to Figure 8(a) obtained by using the semi-global stereo matching algorithm is shown in Figure 8(b), and the disparity map corresponding to Figure 8(a) obtained by using the method of the present invention is shown in Figure 8(c). It can be seen that the disparity optimization effect of the present invention is better.

[0226] The above calculation examples of the present invention are only used to explain the calculation model and calculation process of the present invention in detail, and are not intended to limit the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still within the scope of protection of the present invention.

Claims

1. An unsupervised optimization method for building surface parallax based on local gradient descent, characterized in that: The method specifically comprises the following steps: Step 1: collect the left view and right view of the urban building from different perspectives, obtain a disparity map based on the left view and the right view, and perform superpixel segmentation on the pixel points in the disparity map based on the information in the satellite image of the urban building and the disparity map to obtain a superpixel segmentation result; Step 2: Initialize the number of iterations q = 1; Step 3, initialize the number of superpixels k=1; Step 4: Calculate the parameters of the homography model H according to the disparity map corresponding to the kth superpixel; Step 5: performing bilateral filtering on the disparity map corresponding to the k-th superpixel to obtain a disparity map after bilateral filtering; Step 6: According to the disparity map after bilateral filtering, the matching point (x′, y′) corresponding to the pixel point (x, y) in the left view in the right view is re-determined, (x, y) and (x′, y′) are the re-determined matching point pairs, and the loss function corresponding to the pixel point (x, y) is calculated according to the homography model parameters and the re-determined matching point pairs; Similarly, the loss function value corresponding to each pixel in the left view is obtained respectively; Step 7: Determine whether k=K is satisfied; If satisfied, proceed to step eight; Otherwise, set k=k+1 and return to step 4; Step 8: Determine whether the loss function values ​​corresponding to each pixel of the disparity map are less than the disparity threshold ε disp ; If the loss function values ​​corresponding to each pixel in the disparity map are all less than the disparity threshold, the disparity map obtained in the previous iteration is the final disparity map; If the loss function values ​​corresponding to the pixels in the disparity map are not all less than the disparity threshold, execute step nine; Step nine, determine whether the number of iterations q is equal to the threshold Q; If the number of iterations q is equal to the threshold Q, the x′ coordinate of the matching point (x′, y′) is updated according to the loss function value corresponding to each pixel point, and the final disparity map is obtained according to the updated coordinates; If the number of iterations q is not equal to the threshold Q, the x′ coordinate of the matching point (x′, y′) is updated according to the loss function value corresponding to each pixel point, and the disparity map is updated according to the updated coordinates, and then q=q+1 is set, and the process returns to step 3.

2. The unsupervised optimization method for building surface parallax based on local gradient descent according to claim 1, characterized in that: In the step 1, super-pixel segmentation is performed on the pixels in the disparity map according to the information in the satellite image of the urban building and the disparity map to obtain a super-pixel segmentation result; the specific process is: Step 1: Initialize the number of superpixels to K; Initialize the label matrix Met label The size of the label matrix is ​​the same as that of the disparity map. Each element value in the label matrix represents the superpixel to which the pixel at the corresponding position in the disparity map belongs, and the label matrix Met label The initial value of each element in is -1; Initialize the distance matrix Met distance The same size as the disparity map, the distance matrix Met distance Each element value in represents the distance from the pixel point at the corresponding position in the disparity map to the seed point of the superpixel to which it belongs, and the distance matrix Met distance The initial value of each element in is infinite; Then the initial size S of each superpixel is: Where N is the total number of pixels in the disparity map; Step 12: Divide the entire disparity map into K super-pixel regions of size S, and the width and height of each super-pixel region are equal; Step 13: Map each pixel in the disparity map to a 6-dimensional vector. For pixel i, map pixel i to a vector (x i ,y i ,L i ,a i ,b i ,d i ); Among them, x i ,y i represents the plane space coordinates of pixel i, L i ,a i ,b i represents the color coordinates of pixel i in the satellite image mapped to the CIELab color space, d i Represents the disparity value of pixel i in the disparity map; Step 14: For any superpixel region, calculate the mean of the vectors corresponding to all pixels in the superpixel region to obtain the mean vector, and use the mean vector as the vector corresponding to the initial seed point of the superpixel region (x cen ,y cen ,L cen ,a cen ,b cen ,d cen ); in, S 2 is the total number of pixels in the superpixel area; Similarly, the vectors corresponding to the initial seed points of each superpixel region are obtained respectively; Step 15, initialize the number of iterations l = 1; Step 16: for pixel point i in the disparity map, determine whether pixel point i exists in the 2S×2S area around each seed point, determine all seed points in the surrounding 2S×2S area where pixel point i exists, and then calculate the comprehensive distance between the vector corresponding to pixel point i and the vector corresponding to each determined seed point, and obtain the minimum comprehensive distance; If the minimum comprehensive distance obtained is less than the pixel i in the distance matrix Met distance (x, y), then use the obtained minimum comprehensive distance to update the pixel i in the distance matrix Met distance The corresponding distance in (x, y) is calculated, and the pixel i is assigned to the superpixel corresponding to the minimum comprehensive distance. The label matrix Met is updated according to the superpixel to which the pixel i belongs. label The pixel value corresponding to the pixel i in ; Similarly, each pixel in the disparity map is processed separately; Step 17: For any superpixel, calculate the mean of the vectors corresponding to all the pixels belonging to the superpixel after the current iteration, and use the calculated mean vector as the vector corresponding to the new seed point of the superpixel; Similarly, the vector corresponding to the new seed point of each superpixel is calculated respectively; Step 18: Determine whether the number of iterations l is greater than 20, and for each superpixel, calculate the comprehensive distance between the vector corresponding to the new seed point obtained in the current iteration and the vector corresponding to the seed point obtained in the previous iteration; If the number of iterations is greater than 20 or the comprehensive distances corresponding to the seed points of all superpixels are less than the threshold ε sup , then the final superpixel segmentation result is obtained and the final label matrix Met is returned label ; Otherwise, set the number of iterations l=l+1, and use the vectors corresponding to the new seed points of each superpixel to return to step 16.

3. The unsupervised optimization method for building surface parallax based on local gradient descent according to claim 2, characterized in that: The comprehensive distance between the vector corresponding to the pixel point i and the vector corresponding to each determined seed point is calculated respectively, specifically: in, is the color space distance, is the plane space distance, is the parallax space distance; Through and Perform weighted summation to calculate the comprehensive distance D i : Among them, α, β and γ represent and , and α+β+γ=1.

4. The unsupervised optimization method for building surface parallax based on local gradient descent according to claim 1, characterized in that: The specific process of step 4 is as follows: Step 41: Use the homography model H to characterize the building surface model. The homography model H is specifically: Among them, h i′j′ Represents the parameters in the homography model H; Step 4.2: Use the homography model H to transform the matching points (x i ,y i ) and (x i ′,y i ′) is expressed as: After finishing, we get: Rearrange the equation system of formula (8) to obtain the linear equation system: Then, the objective function H(θ) is established according to the linear equation system of formula (9): H(θ)=argMin{err1,err2,...,err N } (10) Where n = 1, 2, ..., N, N is the number of randomly selected matching point pairs, each group includes 4 matching point pairs, and the 4 matching point pairs are non-collinear, N = 4S 2 , the intermediate variable err n The calculation based on the randomly selected nth group of matching points is as follows: The parameters of the homography model are obtained by solving the objective function.

5. The unsupervised optimization method for building surface parallax based on local gradient descent according to claim 1, characterized in that: The specific process of the bilateral filtering process is as follows: Among them, BF[d i ] represents the disparity value d of pixel i in the disparity map i The result of bilateral filtering is represents the neighborhood centered on pixel i, d i and d j Represent the disparity values ​​of pixel i and pixel j respectively, represents the weight of the spatial domain, represents the weight of the disparity domain, represents the normalized weight, ||·|| represents the 2-norm, and |·| represents the absolute value.

6. The unsupervised optimization method for building surface parallax based on local gradient descent according to claim 5, characterized in that: The normalized weight for:

7. The unsupervised optimization method for building surface parallax based on local gradient descent according to claim 5, characterized in that: The spatial domain weight The calculation method is as follows: Among them, (x, y) is the coordinate of pixel i, (u, v) is the coordinate of pixel j, σ s Represents the standard deviation of the spatial domain.

8. The unsupervised optimization method for building surface parallax based on local gradient descent according to claim 5, characterized in that: The disparity domain weight The calculation method is as follows: Among them, σ r Represents the standard deviation of the disparity domain.

9. The unsupervised optimization method for building surface parallax based on local gradient descent according to claim 1, characterized in that: The loss function corresponding to the pixel point (x, y) is calculated based on the parameters of the homography model and the re-determined matching point pair, specifically: Among them, LF(x,y) is the loss function corresponding to the pixel point (x,y).

10. The unsupervised optimization method for building surface parallax based on local gradient descent according to claim 9, characterized in that: The x′ coordinate of the matching point (x′, y′) is updated according to the loss function value corresponding to each pixel point, and the disparity map is updated according to the updated coordinate; the specific process is: Calculate the gradient of the loss function Then the x′ coordinate is updated as follows: Among them, x″ represents the update result of x′, and α represents the learning rate; Then in the updated disparity map, the disparity corresponding to the pixel point (x, y) is: Among them, disp v ' al Indicates the disparity corresponding to the updated pixel point (x, y), disp val =x′-x.

Images