Three-dimensional reconstruction method integrating laser radar and oblique photography

Abstract

The present invention provides a three-dimensional reconstruction method integrating a laser radar and oblique photography, comprising: S1, acquiring an oblique image and image data and preprocessing same to obtain a first oblique image, wherein the image data comprises POS data; S2, acquiring coordinates of a ground control point and high-precision POS data; S3, correcting the oblique image to obtain a second oblique image; S4, performing aerotriangulation on the second oblique image, and obtaining an effective matching point pair by means of improved SIFT; S5, performing combined adjustment on the basis of the effective matching point pair and performing dense matching to obtain a dense point cloud; S6, acquiring a laser point cloud and preprocessing same; S7, unifying the coordinates of the laser point cloud and the dense point cloud, and using improved ICP to obtain a fused point cloud; and S8, constructing a triangulated irregular network and performing texture mapping to generate a three-dimensional model. In the present invention, on the basis of the advantages of high resolution and strong anti-interference capability of a laser radar, defects in an oblique photography three-dimensional model are mitigated, and meanwhile, an oblique image is corrected and an algorithm used in feature matching is improved, thereby improving the precision, integrity and efficiency of the three-dimensional model.

Description

A 3D reconstruction method integrating lidar and oblique photography Technical Field

[0001] The present invention belongs to the field of three-dimensional reconstruction, and in particular relates to a three-dimensional reconstruction method integrating laser radar and oblique photography. Background Art

[0002] Oblique photography is a high-tech technology developed over the past decade in the global photogrammetry field. It transcends the previous limitation of orthophotography, which was limited to vertical imagery. By deploying multiple sensors on a single flight platform, it simultaneously captures images from five different angles—one vertical and four oblique angles—presenting an intuitive, real-world image that matches human perception. It not only accurately reflects the terrain and captures object texture information with high precision, but also generates realistic 3D city models through advanced positioning, fusion, and modeling techniques. This technology has been widely adopted globally in industries such as emergency rescue, urban planning, and surveying and mapping.

[0003] In recent years, 3D reconstruction technology has been widely used in urban planning, land resource surveys, and cultural heritage preservation and restoration. Its application to power grid equipment has also garnered significant attention. For example, 3D reconstruction of transmission towers facilitates the safe operation of power grids and improves the efficiency of transmission line maintenance and inspection. High-precision and efficient acquisition of 3D data has become a research hotspot in real-world 3D modeling. Three primary methods for obtaining 3D reconstruction data sources are manual modeling, drone-based oblique photography, and lidar modeling.

[0004] Manual modeling techniques, with their long production cycles and heavy measurement workload, are suitable for 3D modeling of smaller objects. Otherwise, modeling large-scale transmission towers is neither economical nor effective. 3D reconstruction based on oblique photography, with its high efficiency and low cost, is widely used in various areas of 3D modeling.

[0005] However, oblique photography technology also has certain limitations. Oblique images have large differences in scale, viewing angle, and a large number of blocked areas, which makes dense matching of oblique images very difficult. In addition, due to the limitations of oblique photography shooting angles and operating methods, the three-dimensional point cloud encrypted by aerial triangulation has holes in blocked areas and cannot establish a complete three-dimensional model, which will cause defects in the aesthetics and integrity of the model. Laser point clouds have the advantages of high precision and high density, and a simple laser radar only obtains point cloud information, and the reconstruction obtains the model structure, which cannot truly restore the real scene. In addition,

[0006] Therefore, it is particularly important to design a 3D reconstruction method that integrates lidar and oblique photography and improve the accuracy of 3D reconstruction by leveraging the advantages of both.

[0007] Summary of the Invention

[0008] In order to solve the above technical problems, the present invention provides a three-dimensional reconstruction method that integrates lidar and oblique photography. It utilizes the high resolution and strong anti-interference ability of lidar to optimize the defects of the oblique photography three-dimensional model and solve the problems of ground distortion and holes. In addition, by correcting the oblique images, the accuracy of dense matching is improved, and the precision and integrity of modeling are improved.

[0009] The present invention provides a three-dimensional reconstruction method integrating laser radar and oblique photography, comprising the following steps:

[0010] S1 acquires an oblique image and performs a first preprocessing to obtain a first oblique image and image data, wherein the image data includes POS data; the POS data includes image posture angle information;

[0011] S2 selects image control points on site, obtains their coordinate information, and converts and refines the plane coordinates and elevation values ​​in the POS data to obtain high-precision POS data;

[0012] S3 calculates the transformation matrix between the first oblique image and its corresponding vertical image based on the high-precision POS data, and corrects the first oblique image using the transformation matrix to obtain a second oblique image;

[0013] S4 uses high-precision POS data and image control point coordinate information to perform aerial triangulation on the second oblique image, and uses the improved SIFT algorithm to extract and match feature points of the second oblique image to obtain valid matching point pairs;

[0014] S5 performs multi-view joint adjustment based on valid matching point pairs to obtain a sparse point cloud, and then uses a multi-view dense matching algorithm to generate a dense point cloud.

[0015] S6 obtains the laser point cloud of the laser radar and performs a second preprocessing;

[0016] S7 performs coordinate uniformity conversion on the pre-processed laser point cloud and the dense point cloud, and uses the improved ICP algorithm to perform point cloud registration and fusion to obtain a fused point cloud;

[0017] S8 constructs a triangulated network based on the fused point cloud, smoothes and simplifies the triangulated network, and automatically maps textures based on the oblique photography image to generate a three-dimensional model. Specifically, the correction process in step S3 includes:

[0018] S31 calculates a first rotation matrix of the first tilted image using the image posture angle information in the high-precision POS data;

[0019] S32 simulates the transformation of the image shooting angle of view, and calculates the transformation matrix between the first oblique image and its corresponding vertical image at each angle of view;

[0020] S33 resamples the first oblique image through a transformation matrix, transforms the oblique images at different oblique viewing angles to a vertical viewing angle, and obtains a second oblique image.

[0021] Specifically, step S4 further includes:

[0022] S41 uses Gaussian pyramid to construct the scale space of the second oblique image based on SIFT algorithm;

[0023] S42 uses the MSERs algorithm to perform extreme value detection on the constructed scale space and locates the feature points in the space based on the detection results;

[0024] S43 uses the binary FREAK descriptor to describe the feature points and obtain the feature vector;

[0025] S44 performs rough matching on the feature vectors according to cosine similarity to obtain primary matching point pairs;

[0026] S45 uses the RANSAC algorithm to eliminate incorrect matching point pairs in the primary matching point pairs and obtain valid matching point pairs.

[0027] Specifically, step S7 further includes:

[0028] S71 integrates the dense point cloud and the laser point cloud to obtain an integrated point cloud, calculates the center of gravity of the integrated point cloud, calculates a second rotation matrix and a translation matrix with the center of gravity as the center, and performs rotation and translation respectively to obtain a dense point cloud and a laser point cloud with unified coordinates;

[0029] S72 performs feature extraction on the dense point cloud and the laser point cloud with unified coordinates, calculates the normal vector and umbrella curvature features of each point cloud, and obtains point cloud features;

[0030] S73 matches and fuses dense point clouds and laser point clouds with the same point cloud features using the ICP algorithm based on the bidirectional KD tree.

[0031] Specifically, the multi-view dense matching algorithm in step S5 is a multi-view dense matching algorithm based on clustering and patch models.

[0032] Specifically, the first preprocessing in step S1 includes light and color uniformity and distortion removal.

[0033] Specifically, the plane coordinates in step S2 are transformed using Gauss-Krüger rectangular plane projection; the elevation values ​​are transformed using the seven parameters obtained from the image control points, or are converted using a quasi-geoid refined surface.

[0034] Specifically, the second preprocessing in step S6 includes denoising and smoothing.

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

[0036] (1) The present invention integrates the high efficiency of drone oblique photography and the accuracy of laser radar, solves the problems of ground distortion and holes in modeling with a single data source, and improves the accuracy, integrity and efficiency of modeling; and adopts a bidirectional KD tree method to match dense point clouds and laser point clouds, thereby improving the accuracy of feature point extraction, providing accurate feature matching data for point cloud fusion, and thus improving the accuracy of three-dimensional reconstruction.

[0037] (2) The present invention calculates the transformation matrix between the oblique image and its corresponding vertical image through the correspondence between the orientation elements and the oblique image, corrects the oblique image through the transformation matrix, transforms the image at different oblique angles to the vertical angle, improves the accuracy of feature extraction and matching of the oblique image, and thus improves the accuracy and integrity of modeling;

[0038] (3) The present invention addresses the problems of the SIFT algorithm, such as long feature point detection time, large amount of calculation, and low precision and poor effect of tilt image matching. Based on the SIFT algorithm, the present invention adopts the MSERs algorithm as the feature point detection algorithm and the FREAK algorithm as the feature point descriptor, adopts a coarse-to-fine matching strategy, improves the precision and efficiency of image matching, and thus improves the efficiency and precision of 3D reconstruction.

[0039] (4) The present invention combines the efficient algorithm of the KD tree search neighborhood with the ICP algorithm to accurately align and fuse point clouds, greatly improving the efficiency of aligning massive point clouds, thereby further improving the efficiency and accuracy of three-dimensional reconstruction. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0041] FIG1 is a schematic diagram of a method according to an embodiment of the present invention;

[0042] FIG2 is a flowchart of a correction process for a first tilted image according to an embodiment of the present invention;

[0043] FIG3 is a technical flow chart of generating a three-dimensional model according to an embodiment of the present invention; Specific embodiments

[0044] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention are within the scope of protection of the present invention.

[0045] Please refer to FIG1 , which is a schematic diagram of a method according to an embodiment of the present invention;

[0046] The present invention provides a three-dimensional reconstruction method integrating laser radar and oblique photography, comprising the following steps:

[0047] S1 acquires an oblique image and performs a first preprocessing to obtain a first oblique image and image data, wherein the image data includes POS data; the POS data includes image posture angle information;

[0048] A five-lens camera is mounted on the same flight platform, capturing ground images from both vertical and oblique angles. The oblique images include a set of positive images taken with the camera perpendicular to the ground, and four sets of oblique images taken with the camera facing at a certain angle to the ground.

[0049] Oblique photography technology is different from traditional vertical aerial photography not only in terms of photography methods, but also in its post-processing data and results. The main purpose of oblique photography technology is to obtain information on multiple directions (especially the side) of the ground objects and to allow users to browse from multiple angles, conduct real-time measurements, browse in three dimensions, and obtain various information. The image data carries the position information (GPS) of the drone at the time of shooting, that is, the line elements x, y, z and attitude information at the moment of camera exposure, and the flight attitude information in the positioning and attitude determination system (POS), that is, the attitude angle elements recorded at the moment of camera exposure. ω, k.

[0050] Oblique photography flight paths are designed using dedicated flight path planning software, ensuring that relative altitude, ground resolution, and physical pixel size maintain trigonometric proportions. Flight paths are typically designed with a 30% lateral overlap and a 66% directional overlap. Currently, for automated model production, both 66% lateral overlap and 66% directional overlap are required. The flight path planning software generates a flight plan file containing the aircraft's flight path coordinates and the exposure point locations for each camera. During flight, each camera automatically exposes and captures images based on the corresponding exposure point coordinates.

[0051] Oblique images have the following characteristics: They reflect the real situation around the objects. Compared with orthophotos, oblique images allow users to observe objects from multiple angles, reflecting the actual situation of the objects more realistically, which greatly makes up for the shortcomings of applications based on orthophotos. Oblique images can realize single image measurement. Through the application of supporting software, measurements including height, length, area, angle, slope, etc. can be directly performed based on the result image, which expands the application of oblique photography technology in the industry. The texture of the side of the building can be collected. For various three-dimensional digital city applications, the large-scale mapping characteristics of aerial photography are utilized, and the batch extraction and texture pasting methods from oblique images can effectively reduce the cost of three-dimensional modeling. It is easy to publish on the Internet. The data format of oblique images can be quickly published on the Internet using mature technology to achieve shared applications.

[0052] In an embodiment of the present invention, the first preprocessing in step S1 includes light and color uniformity and distortion removal.

[0053] Specifically, due to temporal and spatial differences during flight, color shifts can occur between images, necessitating uniform light and color processing. The image data acquired from each flight is promptly and carefully inspected and preprocessed, and the aerial images are adjusted and generated strictly according to the uniform light and color steps, ultimately obtaining image data (in TIFF format) with optimal imaging effects. All results undergo self-quality inspections and are archived after passing the inspection, resulting in the final aerial photography results.

[0054] Image preprocessing is a crucial step in transforming aerial images from invisible to visible and achieving color restoration. In digital aerial photography, image preprocessing's impact on post-production results lies in processing speed and adjusting for uniform lighting and color. To obtain the best data, uniform lighting and color processing of the oblique images involves: distributed, uninterrupted data processing to ensure the fastest possible data processing, improve feedback efficiency, and expedite project completion. The principle of image preprocessing is to ensure that the image histogram covers the entire range of 0-255 color levels as closely as possible and approximates a normal distribution. This ensures that true-color images have rich tones, saturated colors, good color balance, and normal color reproduction. In aerial photography conditions with minimal light differences, the same color grading template is used for the same route. In situations with large light differences and low atmospheric transparency, individual color grading is performed on one or more images to achieve optimal, realistic colors. The oblique camera is equipped with five cameras facing in different directions, with a 45° oblique shooting angle. At the moment of exposure, the five cameras capture images from different angles, inevitably leading to inconsistencies in light contrast and intensity. This can cause variations in color and brightness for the same feature imaged by different cameras, ultimately impacting the accuracy and quality of the 3D real-world scene. Drawing on previous experience with oblique photogrammetry data processing, image color can be adjusted during the image quality inspection and mosaic stages to address issues such as fog and poor contrast in parts of the test area caused by weather, eliminating the effects of these factors. Following pre-processing of the aerial photography data, all image data undergoes a further light and color uniformity process to further unify the image color across the entire area. This reduces color grading time in later production and ensures a more natural and uniform color throughout the project. Aerial photography data requires post-processing to provide high-quality aerial image data after light uniformity, color grading, and distortion removal. Data processing records and reports are also provided to meet the standards and specifications for data production.

[0055] Please refer to FIG2 , which is a flowchart of a correction process for a first tilted image according to an embodiment of the present invention;

[0056] S2 selects image control points on site, obtains their coordinate information, and converts and refines the plane coordinates and elevation values ​​in the POS data to obtain high-precision POS data;

[0057] In an embodiment of the present invention, the plane coordinates in step S2 are transformed using Gauss-Krüger rectangular plane projection; the elevation values ​​are transformed using the seven parameters obtained from the image control points, or are converted using a quasi-geoid refined surface.

[0058] S3 calculates the transformation matrix between the first oblique image and its corresponding vertical image based on the high-precision POS data, and corrects the first oblique image using the transformation matrix to obtain a second oblique image;

[0059] In an embodiment of the present invention, the correction process in step S3 includes:

[0060] S31 calculates a first rotation matrix of the first tilted image using the image posture angle information in the high-precision POS data;

[0061] S32 simulates the transformation of the image shooting angle of view, and calculates the transformation matrix between the first oblique image and its corresponding vertical image at each angle of view;

[0062] S33 resamples the first oblique image through a transformation matrix, transforms the oblique images at different oblique viewing angles to a vertical viewing angle, and obtains a second oblique image.

[0063] In the embodiment of the present invention, the image posture angle information in the high-precision POS data (taking the OPK angle as an example) is brought into the following formula to calculate the rotation matrix R N ;

[0064] in, a1, a2, a3, b1, b2, b3, c1, c2 and c3 are the 9 direction cosines at the moment of camera exposure; the direction cosines at the moment of camera exposure and the image attitude angle elements ω, k;

[0065] According to R i =R N R Ni Calculate the first rotation matrix R of the oblique image i , and at the same time by t i =R N t Ni +t N Calculate the line element t of the oblique image i , where i and N represent oblique image and vertical image respectively.

[0066] According to the relationship between image point and object point in projection transformation, we can get: δ·X=K·(R cv ·(xyz) T +T)

[0067] Where δ is the scale factor, K is the camera matrix, T=-R cv *t, X=(xyz) T , where (xyz) is the three-dimensional coordinate of the image point.

[0068] Let z=0, the above formula becomes: δ·X=K·(R cv ·(xy 0) T + T) = K·([r1 r2 r3]·(xy 0) T +T) =K·[r1 r2 T]·(xy 1) T

[0069] Let R kt =K·[r1 r2 T], then δ·X=R kt ·(xy 1) T .

[0070] Then the corresponding relationship between the vertical image point and the object point is as follows: δ N ·X N =R Nkt ·(xy 1) T

[0071] Similarly, the corresponding relationship between the image point and the object point of the oblique viewing angle image is as follows: δ i ·X i =R ikt ·(xy 1) T

[0072] Combining the above relations, we get:

[0073] In the above formula, the matrix R cst That is, the transformation matrix corresponding to the oblique and vertical images is used to correct the images of each perspective in the oblique image to obtain the corrected oblique image.

[0074] Please refer to FIG3 , which is a technical flow chart of generating a three-dimensional model according to an embodiment of the present invention;

[0075] S4 uses high-precision POS data and image control point coordinate information to perform aerial triangulation on the second oblique image, and uses the improved SIFT algorithm to extract and match feature points of the second oblique image to obtain valid matching point pairs;

[0076] In an embodiment of the present invention, step S4 is specifically as follows:

[0077] S41 uses Gaussian pyramid to construct the scale space of the second oblique image based on SIFT algorithm;

[0078] Scale-Invariant Feature Transform (SIFT) is an algorithm used to detect and describe local features of images, find extreme points in spatial scale, and extract their position, scale, and rotation invariants.

[0079] The main steps of extracting feature points using the improved SIFT algorithm are as follows:

[0080] 1) Constructing Gaussian scale space

[0081] In order to find stable feature points in different scale spaces, the SIFT algorithm uses Gaussian difference kernels of different scales to generate Gaussian difference scale space DOG. The specific formula is as follows: D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y)=L(x,y,kσ)-L(x,y,σ)

[0082] The Gaussian convolution kernel is the only linear kernel that can achieve scale transformation, so the scale space of a two-dimensional image is defined as: L(x,y,σ)=G(x,y,σ)*I(x,y)

[0083] Where (x, y) is the spatial coordinate of the image pixel, I(x, y) represents the pixel value of the original image, the size of σ determines the smoothness of the image, G(x, y, σ) is a scale-variable Gaussian function, k is the scale space factor, and L(x, y, kσ) is the Gaussian function at the corresponding scale.

[0084] S42 uses the MSERs algorithm to perform extreme value detection on the constructed scale space and locates the feature points in the space based on the detection results;

[0085] The MSERs algorithm is similar to the watershed concept. It first extracts the grayscale image of the original image. Generally, the grayscale range of the grayscale image is 0 to 255. Then, a threshold is set within this range. The threshold increases from 0 to 255. The grayscale value on the two-dimensional image can be regarded as its "height". When the threshold increases, the pixels with low grayscale in the image are "flooded", and the grayscale value of the image in this area becomes 0. This process is similar to the process of continuously injecting water into a hilly area, which is the watershed concept.

[0086] MSERs have affine-invariant stability and are invariant to image grayscale changes with affine transformations. The support set of this region is stable relative to grayscale changes and can be used to detect regions of varying degrees of detail. When the image is binarized using different grayscale thresholds, the most stable region is obtained, and the connected regions between adjacent threshold images are ultimately MSERs. Therefore, MSERs are divided into MSERs+ and MSERs-. For two images with affine transformations obtained from the same scene, the extracted radially invariant regions differ in distortion, scale, and rotation. Although MSERs obtains an affine-invariant region fitted with an ellipse, in order to facilitate subsequent processing, the ellipse must be fitted and normalized to a circle, and the MSER region must be enlarged.

[0087] The specific process is summarized as follows:

[0088] (1) Image grayscale: The input digital image is grayscaled, and then each pixel in the image is sorted from large to small or in reverse order according to the Bin-Sort algorithm;

[0089] (2) Image pixel clustering: Based on the previous step, the sorted image points are mapped to the operation domain through Union-Find, and then the position of each pixel mapping is recorded. The mapping area is searched from bottom to top, and the image pixels are binarized according to the threshold value that is constantly increasing. The pixel area merged during the processing is recorded as Q, where Q represents any pixel point. When p∈Q, q∈σQ, and I(p)>I(q) is satisfied, Q is judged as the maximum extreme value area, otherwise it is the minimum extreme value area, where σQ represents the area boundary of Q.

[0090] (3) Maximum stable extreme value region detection: Based on the extreme value regions obtained in the previous step, such as Q1, ...Q i-1 , Qi, etc., when the above areas meet When the conditions are met, the extreme area change rate q(i)=|Q i+Δ / Q i-Δ | / |Q i |, where |·| is used to calculate the area of ​​the extreme value region, Δ∈{0,1,...,225}, if the result of the above formula calculated by i* is the smallest, then it can be determined This extreme value region is the maximum stable region. Calculate the value of q(i) for each extreme value region that meets the conditions until all the extreme value regions that meet the conditions are The detected area is recorded as MSERs+; then the formula I inver =I max -I performs grayscale inversion processing on the input original image, and then repeats the above two steps, and the extracted extreme value regions are recorded as MSERs-;

[0091] (4) Output of MSERs ellipse region result: The MSERs region extracted by the above steps is irregular. In order to facilitate the description of features, an ellipse is usually used to fit the affine invariant region. The region fitted by the ellipse includes the major and minor semi-axes of the ellipse, i.e., a and b, the inclination angle of the ellipse, and the center of gravity (x g ,y g ). The steps of ellipse fitting are as follows:

[0092] In the simulated region Q, μ pq =ΣΣx p x qI(x,y) defines its (p+q)-order two-dimensional geometric moment, where p,q=0,1,2.... μ 00 =ΣΣI(x,y) μ 10 =∑∑xI(x,y) μ 01 =∑∑yI(x,y) μ 20 =∑∑(xx g ) 2 I(x,y) μ 02 =∑∑(yy g ) 2 I(x,y) μ 11 =∑∑(xx g ) 2 (yy g )I(x,y)

[0093] Among them, x g =μ 10 / μ 00 The calculation formula for the horizontal coordinates of the ellipse is represented by g =μ 01 / μ 00 It is the formula for calculating the vertical coordinate.

[0094] According to the principle steps of MSERs, it is easy to know that the grayscale values ​​of pixels in the extreme value area Q are normalized to 1, so the grayscale values ​​of pixels I(x,y) in the above formula are all 1. Therefore, the covariance of the matrix is The matrix Ω can uniquely determine a characteristic ellipse, and the major and minor semi-axes of the ellipse and their orientation angles can be calculated using the following formula:

[0095] Among them, γ1 and γ2 are the two eigenvalues ​​of the covariance matrix Ω, which can be obtained as follows:

[0096] By taking the above steps on the extracted irregular extreme value region, the fitted MSERs elliptical region can be obtained, thereby detecting all the feature points in the region, which is convenient for the subsequent feature region description.

[0097] S43 uses the binary FREAK descriptor to describe the feature points and obtain the feature vector;

[0098] In an embodiment of the present invention, the feature points determined in step S42 are encoded, and the feature points are described using the FREAK descriptor to obtain feature vectors. The specific method is as follows:

[0099] With the feature point obtained in step S42 as the center, the sampling points are evenly distributed on 7 concentric circles with the feature point as the center. The smaller the interval with the feature point, the denser the sampling points; the larger the interval with the feature point, the sparser the sampling points. In addition, each sampling point needs to be Gaussian smoothed to eliminate the influence of noise on the sampling point. According to the intensity relationship between the sampling point pairs, the binary feature vector is obtained by encoding, as follows: F = ∑0≤a<N2 a T(P a )

[0100] Among them, P a Represents a pair of sampling points, N is the length of the descriptor, T(P a ) is the output result of the sampling point pair.

[0101] in: is the grayscale value intensity of the sampling point after Gaussian smoothing, r1 and r2 are the radii of the concentric circles where different sampling points are located, and each feature point itself participates in the comparison, and finally a feature point descriptor of dimension 2 is obtained.

[0102] Because the sampling points may overlap in some areas during selection, resulting in low discriminability in the original feature vector, the original feature vector needs to be reconstructed, retaining dimensions with greater information content and reducing feature vector redundancy. The iris feature vector is reconstructed using the following steps, selecting and retaining dimensions with high information content: First, a matrix is ​​created for the detected feature points. Each row in the matrix represents a feature point, and each column in each row represents the comparison results of all pairs of sampling points including that feature point. The matrix size is N × 903, where N is the number of detected feature points. Next, the mean of each column in the matrix is ​​calculated. To obtain highly discriminative features, the variance of the column vector should be large, meaning the mean of the corresponding column should be closer to 0.5. Finally, the difference between the mean of each column and 0.5 is calculated. The first 512 columns are sorted by the difference from smallest to largest, resulting in the feature vector.

[0103] S44 performs rough matching on the feature vectors according to cosine similarity to obtain primary matching point pairs;

[0104] In the embodiment of the present invention, cosine similarity, i.e., cosine distance, uses the cosine value of the angle between two feature vectors to calculate their similarity, thereby measuring the size of the difference between two individuals. Compared with Euclidean distance, it pays more attention to the difference in direction between the two vectors. Therefore, the present invention uses cosine similarity to perform rough matching of feature points, assuming that the n-dimensional feature vectors of the primary matching point pair are A and B respectively; A is [A1, A2, ..., A n ], where B is [B1, B2, ..., B n]; AT represents the transposed matrix of matrix A, and the cosine similarity of the primary matching point pair is represented by cosθ:

[0105] By calculating the cosine value, i.e., the value of cosθ, between the feature vectors of the above matching point pairs, when the value of cosθ is ≥ 0.9, it is considered that the similarity of the matching point pair in direction is also relatively high, and it is considered to be a correct matching point pair. Otherwise, it is an incorrect matching point pair and is discarded.

[0106] S45 uses the RANSAC algorithm to eliminate incorrect matching point pairs in the primary matching point pairs and obtain valid matching point pairs.

[0107] In an embodiment of the present invention, the purpose of eliminating mismatched point pairs is achieved by using the Random Sample Consensus Algorithm (RANSAC) supplemented by dynamically adjusting the ratio of the nearest neighbor to the next nearest neighbor.

[0108] S5 performs multi-view joint adjustment based on valid matching point pairs to obtain sparse point cloud and corrected camera parameters; the sparse point cloud is converted into a dense point cloud using a multi-view dense matching algorithm.

[0109] In an embodiment of the present invention, the multi-view dense matching algorithm in step S5 is a multi-view dense matching algorithm based on clustering and patch models.

[0110] Based on the sparse point cloud, the Clustering Multiple View Stereo & Patch-based Multi-View Stereo (CMVS+PMVS) algorithm is used to generate a dense point cloud. That is, based on the sparse point cloud, the patches are divided and the patch structure is expanded through patch diffusion to obtain a dense point cloud.

[0111] S6 obtains the laser point cloud of the laser radar of the UAV flight target area and performs a second preprocessing;

[0112] In an embodiment of the present invention, the second pre-processing in step S6 includes denoising and smoothing, wherein the laser point cloud is subjected to denoising and smoothing to filter out some obvious outliers.

[0113] S7 performs coordinate uniformity conversion on the pre-processed laser point cloud and the dense point cloud, and uses the improved ICP algorithm to perform point cloud registration and fusion to obtain a fused point cloud;

[0114] In an embodiment of the present invention, step S7 is specifically as follows:

[0115] S71 integrates the dense point cloud and the laser point cloud to obtain an integrated point cloud, calculates the center of gravity of the integrated point cloud, calculates a second rotation matrix and a translation matrix with the center of gravity as the center, and performs rotation and translation respectively to obtain a dense point cloud and a laser point cloud with unified coordinates;

[0116] The dense point cloud and the laser point cloud are integrated to obtain the integrated point cloud. The center of gravity of the integrated point cloud is calculated. All point clouds are rotated and translated around the center of gravity. The rotation angle is the angle at which each point cloud is rotated relative to the initial point cloud.

[0117] The prerequisite for data fusion is to ensure the consistency of the dense data space coordinate systems of the laser point cloud and oblique image. The following two calculation formulas can be used to solve the rotation and translation parameters required for the coordinate transformation of the two point cloud data.

[0118] Where (x, y, z) is the initial coordinates of each point in the integrated point cloud, and (X, Y, Z) is the coordinates of the point cloud after coordinate transformation. α, β, and γ are rotation parameters, and Δx, Δy, and Δz are translation parameters. By manually selecting the 3D coordinates of feature points with the same name, the second rotation matrix and translation matrix can be calculated, resulting in the six parameters α, β, γ, Δx, Δy, and Δz.

[0119] S72 performs feature extraction on the dense point cloud and the laser point cloud with unified coordinates, calculates the normal vector and umbrella curvature features of each point cloud, and obtains point cloud features;

[0120] First, the least squares method is used to estimate the normal vector of each point cloud data through plane fitting; then the umbrella curvature of each point cloud data is calculated: let P0 be the center point, n is the normal vector of P0 point, N i (i=1,2,…,8) are the neighborhood points of P0, then the curvature of the center point P0 is:

[0121] Among them, ABS means taking the absolute value, n i is the neighborhood point N i The normal vector of

[0122] Imagine an umbrella in three states: fully open, open, and closed. When the umbrella is closed, the projected area of ​​the umbrella, calculated by combining the vector product of the center point and the eight neighboring points, is at its maximum, approaching 8. As the umbrella opens, the value decreases until the fully open projected area approaches zero, and the umbrella's curvature also approaches zero.

[0123] S73 matches and fuses dense point clouds and laser point clouds with the same point cloud features based on the ICP algorithm of the bidirectional KD tree to obtain a fused point cloud.

[0124] During the point cloud fusion process, an ICP algorithm based on a bidirectional KD tree is used, taking into account the point cloud normal vector and curvature. When calculating the target error, if the point pairs are all feature points, they are assigned a larger weight. This method is used to select more accurate point pairs for point cloud fusion. The present invention uses an improved ICP algorithm for precise registration of image-based point clouds. When fusing dense point clouds and laser point clouds, a KD tree spatial search method is used to fully utilize the relationship between data points and their neighboring points, efficiently process the point cloud data, and provide accurate feature matching points for point cloud registration.

[0125] The KD tree algorithm is a multidimensional search tree structure that extends the binary tree search to three-dimensional space. The KD tree construction algorithm uses a step-by-step recursive method to determine the splitting dimension based on the maximum variance. The node where the splitting dimension is located is used as the root node to divide the data structure into two left and right spaces. The left and right spaces are then recursively divided until each subset no longer contains point clouds. This method greatly improves the efficiency of precise point cloud registration, thereby further improving the efficiency of three-dimensional reconstruction.

[0126] In an embodiment of the present invention, an ICP algorithm based on a bidirectional KD tree is used to achieve accurate registration of point clouds. A bidirectional KD tree method is used to search for neighboring domain points in cloud data. The matching and fusion steps are as follows:

[0127] S731 For the dense point cloud R and the laser point cloud S, use the KD tree method to find the matching points of P in the point cloud S to form the first initial corresponding point set (r, s), where r corresponds to the point in the point cloud R and s corresponds to the point in the point cloud S;

[0128] S732 finds matching points in point cloud S and its corresponding points in point cloud R to form a second initial corresponding point set (r', s'), where r' corresponds to a point in point cloud P and s' corresponds to a point in point cloud S;

[0129] S733 takes the intersection of the first initial corresponding point set (r, s) and the second initial corresponding point set (r', s') to obtain valid feature points of the two sets of point clouds, and starts iteration with the valid feature points as the final corresponding point set for point cloud matching to obtain the matching relationship between the final dense point cloud and the laser point cloud, and fuses them using the matching relationship to obtain a fused point cloud.

[0130] The intersection of the initial corresponding points is considered the exact correspondence between the two point clouds. This relationship is used to adjust the overall position and angle of the two point clouds, and finally form a single point cloud. The result of fusion is that the two point clouds become one point cloud, and finally the point clouds from multiple angles are merged into one, which means the number of points increases several times.

[0131] The bidirectional KD tree method can effectively extract point cloud features and improve the accuracy of point cloud registration. The ICP algorithm of the bidirectional KD tree is used for cyclic iteration, taking the point cloud normal vector and curvature into account. When calculating the target error, if all point pairs are feature points, a larger weight is set for them. This method is used to select more accurate point pairs for point cloud fusion to obtain the final matching relationship between the dense point cloud and the laser point cloud.

[0132] S8 builds a triangulated network based on the fused point cloud, smoothes and simplifies the triangulated network, and automatically maps textures based on oblique photography images to generate a three-dimensional model.

[0133] The process for constructing an irregular triangulated network is as follows: A terrain dataset is converted to an irregular triangulated network using the "Terrain to TIN" geoprocessing tool. Texture mapping involves mapping texture pixels from texture space to screen space. First, a texture object is defined. Then, the glBlindTexture function is used to select the texture object and load the corresponding texture when rendering the 3D scene. Finally, the glBlindTexture function is used to delete the texture object before outputting the result.

[0134] The present invention adopts three-dimensional lidar technology to assist oblique photogrammetry to repair the missing point clouds in blind areas, so that the two can give full play to their respective advantages and characteristics, ensuring both production efficiency and the integrity of the point clouds in the obscured areas of oblique photography, and ultimately improving the fineness of the model.

[0135] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present invention should be included in the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be based on the scope of protection of the claims.

Claims

1. A three-dimensional reconstruction method integrating laser radar and oblique photography, characterized in that: The following steps are involved: S1 acquires an oblique image and performs a first preprocessing to obtain a first oblique image and image data, wherein the image data includes POS data; the POS data includes image posture angle information; S2 selects image control points on site, obtains the coordinate information of the image control points, and converts and refines the plane coordinates and elevation values ​​in the POS data to obtain high-precision POS data; S3 calculates the transformation matrix between the first oblique image and its corresponding vertical image according to the high-precision POS data, and corrects the first oblique image by the transformation matrix to obtain a second oblique image; S4 performs aerial triangulation on the second oblique image through high-precision POS data and image control point coordinate information, extracts and matches feature points of the second oblique image through an improved SIFT algorithm, and obtains effective matching point pairs; S5 performs multi-view joint adjustment based on valid matching point pairs to obtain a sparse point cloud, and then uses a multi-view dense matching algorithm to generate a dense point cloud. S6 obtains the laser point cloud of the laser radar and performs a second preprocessing; S7 converts the pre-processed laser point cloud and dense point cloud into coordinate uniformity, and uses the improved ICP algorithm to perform point cloud registration and fusion to obtain a fused point cloud; S8 builds a triangulated network based on the fused point cloud, smoothes and simplifies the triangulated network, and automatically maps textures based on oblique photography images to generate a 3D model; The correction process in step S3 includes: S31 calculates a first rotation matrix of the first tilted image through the image posture angle information in the high-precision POS data; S32 simulates the transformation of the image shooting angle of view, and calculates the transformation matrix between the first oblique image and its corresponding vertical image at each angle of view; S33 resamples the first oblique image by using a transformation matrix, transforms the oblique images at different oblique viewing angles to a vertical viewing angle, and obtains a second oblique image; Wherein, step S4 further comprises: S41 uses Gaussian pyramid to construct the scale space of the second oblique image based on SIFT algorithm; S42 uses the MSERs algorithm to perform extreme value detection on the constructed scale space and locates the feature points in the space according to the detection results; S43 uses the binary FREAK descriptor to describe the feature points and obtain the feature vector; S44 performs rough matching on the feature vectors according to cosine similarity to obtain primary matching point pairs; S45 uses the RANSAC algorithm to remove mismatched point pairs in the primary matching point pairs to obtain valid matching point pairs; Wherein, step S7 further comprises: S71 integrates the dense point cloud and the laser point cloud to obtain an integrated point cloud, calculates the center of gravity of the integrated point cloud, calculates a second rotation matrix and a translation matrix with the center of gravity as the center, and performs rotation and translation respectively to obtain a dense point cloud and a laser point cloud with unified coordinates; S72 performs feature extraction on the dense point cloud and the laser point cloud with unified coordinates, calculates the normal vector and umbrella curvature features of each point cloud, and obtains point cloud features; S73 matches and fuses dense point clouds and laser point clouds with the same point cloud features based on the ICP algorithm of the bidirectional KD tree.

2. The three-dimensional reconstruction method integrating laser radar and oblique photography according to claim 1, characterized in that: The multi-view dense matching algorithm in step S5 is a multi-view dense matching algorithm based on clustering and patch models.

3. The three-dimensional reconstruction method integrating laser radar and oblique photography according to claim 1, characterized in that: The first preprocessing in step S1 includes light and color uniformity and distortion removal.

4. The three-dimensional reconstruction method integrating laser radar and oblique photography according to claim 1, characterized in that: The plane coordinates in step S2 are transformed by using Gauss-Krüger rectangular plane projection; the elevation values ​​are transformed by using the seven parameters required by the image control points, or are converted by using a quasi-geoid refined surface.

5. The three-dimensional reconstruction method integrating laser radar and oblique photography according to claim 1, characterized in that: The second preprocessing in step S6 includes denoising and smoothing.

Images