Multi-View Underwater 3D Point Cloud Mosaic Method Based on the Fusion of Sparse and Dense Point Clouds

Abstract

The present invention discloses a multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds, which relates to the technical field of underwater structure detection. The present invention includes: receiving surface images of an underwater structure binocular camera, where the binocular camera includes a left camera and a right camera, generating a dense point cloud based on stereo matching and triangulation, and outputting the correspondence between the dense point cloud and the pixels of the left camera image; based on the received left camera image, using a structure-from-motion algorithm to generate a sparse point cloud by multi-view triangulation, and outputting the correspondence between the sparse point cloud and the pixels of the left camera image. When stitching the point clouds reconstructed from multiple views, the present invention can better restore the surface features of the underwater structure, without the need for feature extraction, has strong accuracy and robustness, and can achieve high-precision point cloud stitching in complex underwater environments. Therefore, it has good application prospects in actual underwater structure detection and image stitching tasks.

Description

Technical Field

[0001] The present invention relates to the technical field of underwater structure detection, and specifically to a multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds. Background Art

[0002] With the wide application of underwater infrastructure (such as bridges, pipelines, offshore platforms, etc.), the structural health monitoring thereof has become the key to engineering safety assurance. However, traditional underwater detection methods mainly rely on manual inspection by divers or single-sensor technology, and have problems such as high operation risks, low detection accuracy, and poor environmental adaptability. In recent years, three-dimensional point cloud technology has gradually received attention in the field of underwater detection. Binocular stereo vision can generate dense point clouds, but it is difficult to obtain stable global features due to underwater lighting and imaging quality. Although the sparse point cloud based on Structure from Motion (SFM) can provide more accurate pose estimation, it lacks sufficient point cloud density. Therefore, the present invention proposes a multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds. Summary of the Invention

[0003] The purpose of the present invention is to provide a multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds, which can achieve high-precision point cloud stitching in complex underwater environments, and the stitching accuracy is close to that in water surface environments, having wide engineering application value.

[0004] According to the first aspect of the present invention, to achieve the above purpose, the present invention provides the following technical solution: A multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds, comprising the following steps:

[0005] Receiving surface images of a binocular camera of an underwater structure, where the binocular camera includes a left camera and a right camera, generating a dense point cloud based on stereo matching and triangulation, and outputting the correspondence between the dense point cloud and the pixels of the left camera image;

[0006] Based on the received left camera image, using a Structure from Motion algorithm to generate a sparse point cloud by multi-view triangulation, and outputting the correspondence between the sparse point cloud and the pixels of the left camera image;

[0007] According to the correspondences between the obtained dense point cloud and sparse point cloud and the image pixels respectively, calculating the correspondence between the dense point cloud and the sparse point cloud, obtaining a preliminary alignment of the point cloud positions, and completing the rough stitching of the point cloud;

[0008] On the basis of the rough stitching of the point cloud, using the Iterative Closest Point algorithm to perform precise stitching on the point cloud, and minimizing the error in the overlapping area of the point cloud.

[0009] Furthermore, the surface image of the underwater structure is acquired by a left camera and a right camera. The two cameras with different directions shoot the surface of the object in the same scene at the same time in different directions, and the morphology of the object surface is reconstructed through the images of the two cameras.

[0010] Furthermore, a dense point cloud is generated based on stereo matching and triangulation, and the correspondence between the dense point cloud and the left camera image pixels is output as follows:

[0011] (31) Use the checkerboard calibration method to calibrate the binocular camera and obtain the intrinsic and extrinsic parameters of the binocular camera:

[0012] (31.1) Camera calibration is the process of establishing the correspondence between the pixel coordinate system and the relative positions of object points in space. It is also the process of determining the camera's intrinsic parameter matrix K, extrinsic parameter matrix Rt, and distortion coefficients. The camera's intrinsic parameter matrix K contains the basic properties of the camera and is used to describe the camera's optical characteristics. Its form is:

[0013]

[0014] in: f x , f y Represents a scaling relationship, through the image distance at the time of imaging f , pixel size s And the actual angle θ of the camera plane pixel arrangement is calculated:

[0015]

[0016] in: s x , s y are the width and height in pixel sizes respectively; (c x , c y ) represents the pixel coordinates of the camera principal point, c x is the pixel x-axis coordinate, c y Pixel y-axis coordinate;

[0017] f s Image distance at imaging f , horizontal pixel size s x And the calculation of the actual angle θ of the sensor plane pixel arrangement:

[0018]

[0019] △r is the mathematical model of radial distortion, expressed as:

[0020]

[0021] The radial distortion is decomposed into the x and y axis directions to obtain a Gaussian radial distortion model. Among them, r is the distance from the image point to the center of the image plane, and r 2 = x 2 + y 2 , k1, k2, and k3 represent the radial distortion coefficients, x and y respectively represent the abscissa and ordinate of the point in the ideal image coordinate system, and △X r represents the Gaussian radial distortion decomposed into the x-axis direction, and △Yr represents the Gaussian radial distortion decomposed into the y-axis direction, which are respectively expressed as:

[0022]

[0023] P(r) is the mathematical model of tangential distortion. Among them, r is the distance from the image point to the center of the image plane, and r 2 = x 2 + y 2 , p1 and p2 represent the tangential distortion coefficients, which are expressed as:

[0024]

[0025] The tangential distortion is decomposed into the x and y axis directions. x and y respectively represent the abscissa and ordinate of the point in the ideal image coordinate system. △x d represents the tangential distortion decomposed into the x-axis direction, and △y d represents the tangential distortion decomposed into the y-axis direction, which are respectively expressed as:

[0026]

[0027] In addition, the external camera parameters include a 3×3 rotation matrix R and a 3×1 translation matrix M to complete the affine transformation from the world coordinate system ( O w - x w y w z w ) to the camera coordinate system ( O - xyz ), where x w , y w , z w respectively represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the point ( x w , y w , z w ) in the three-dimensional world coordinate system.x , y , z respectively represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the point ( x , y , z ) in the three-dimensional camera coordinate system. The affine transformation is expressed as follows:

[0028]

[0029]

[0030] The matrix A is composed of the rotation matrix R and the translation matrix M, where R 11 , R 12 , R 13 , R 21 , R 22 , R 23 , R 31 , R 32 , R 33 respectively represent the corresponding position elements in the rotation matrix R, and M x , M y , M z respectively represent the corresponding position elements in the translation matrix M;

[0031] And in the second step, the camera coordinate system needs to be converted to the image coordinate system ( O i – x i y i ), where x i , y i respectively represent the x-axis coordinate and y-axis coordinate of the point ( x i , y i ) in the two-dimensional image coordinate system. The coordinate conversion is expressed as:

[0032]

[0033] In the formula, f is the image distance during imaging, α is the scale factor of the coordinates, the B matrix is composed of the image distance f . Finally, the image coordinates are converted to the pixel coordinate system ( O s – x s y s ), where x s , ys respectively represent the x-axis coordinate and y-axis coordinate of a point ( x s , y s ) in the pixel coordinate system:

[0034]

[0035] wherein, s x , s y are the width and height of the pixel size; (c x , c y ) represent the pixel coordinates of the camera principal point, θ is the true angle of the pixel arrangement in the camera plane, the C matrix is composed of 1 / s x , -cotθ / s y , c x , 1 / ( s y sinθ), c y . Combining and substituting all coordinate systems, the camera imaging mathematical model is obtained:

[0036]

[0037] ;

[0038] (31.2) The internal parameter matrix K and distortion coefficients of the left camera and the right camera correct the original images, and the two corrected images are in a plane and parallel to each other;

[0039] For the camera coordinate system point (I, J, K), there is:

[0040]

[0041]

[0042]

[0043] where k1, k2, and k3 represent the radial distortion coefficients, p1 and p2 represent the tangential distortion coefficients, and among them, R 11 , R 12 , R 13 , R 21 , R 22 , R 23 , R 31 , R 32 , R 33 respectively represent the corresponding position elements in the rotation matrix R, M x , My , M z They represent the corresponding position elements in the translation matrix M, r is the distance from the image point to the center of the image plane, r 2 =x 2 +y 2 , f x , f y , f s , c x with c y is the element in the internal parameter matrix K, x ’ with y ’ is the horizontal and vertical coordinates of the camera coordinate point after the z coordinate is normalized, x ” with y ” are the horizontal and vertical coordinates of the camera coordinate point after distortion, (x s ,y s ) is the pixel coordinate of the image point, and the intrinsic and extrinsic parameters are calibrated by the checkerboard calibration plate. For each feature point on the calibration plate, the following formula is satisfied:

[0044]

[0045] Where r1, r2, and r3 are the three column vectors of the rotation matrix R, t represents the translation matrix M, K is the internal parameter matrix, and α is the scale factor of the coordinates. w ,y w , z w ) is the coordinate of the point on the chessboard in the world coordinate system, and the calibration point (x w ,y w , z w ) and the image point (x s ,y s ) Direct homography mapping relationship:

[0046]

[0047] Among them, H is the homography matrix, h 11 , h 12 , h 13 , h 21 , h 22 , h 23 , h 31 , h 32 They represent the corresponding position elements in the homography matrix H, s is an unknown scale coefficient, and H ’ =[h 11 h 12 h 13 h 21 h 22 h 23 h31 h 32 T matrix operation equation:

[0048]

[0049] (31.3) Solve the homography matrix H using the least squares method. Since , where h1, h2, and h3 are the three column vectors of the homography matrix H, λ is the scale factor, r1 and r2 are the first two column vectors of the rotation matrix R, and t is the translation matrix T. Since the rotation matrix is an orthogonal matrix, its rows or column vectors are pairwise orthogonal and unit vectors. Thus, two constraint equations are obtained:

[0050]

[0051] K is the camera intrinsic matrix, h1 T , h2 T represent the transposes of the column vectors h1 and h2 respectively. K -T represents the transpose of the inverse matrix of the K vector, and K -1 represents the inverse matrix of the K vector. Let:

[0052]

[0053] B is a symmetric matrix. B 11 , B 12 , B 13 , B 21 , B 22 , B 23 , B 31 , B 32 , B 33 represent the corresponding elements in the matrix B respectively. Define the b matrix with a 6D vector. The T in the upper right corner of the matrix represents the transpose of the matrix:

[0054]

[0055] It can be known that:

[0056]

[0057] h1 and h2 are the first two column vectors of the homography matrix H. h i can take h1 or h2 according to the value of i. h j can take h1 or h2 according to the value of j. h i T represents h i inverse matrix, where v ij T The vector is represented as follows. The T in the upper right corner of the matrix represents the transpose of the matrix:​

[0058]

[0059] h 11 ,h 12 ,h 13 ,h 21 ,h 22 ,h 23 respectively represent the corresponding position elements in the homography matrix H, then the constraint equation is written in matrix operation form using matrix v ij and b:

[0060]

[0061] where v 11 T represents the corresponding position element of v ij T ,v 12 ,v 22 respectively represent the corresponding position elements of v ij obtained from the matrix operation expression. Multiple groups of v ij are substituted into the expression to be written in the form of the solution of the homogeneous equation: Vb = 0; the least squares solution of the homogeneous equation system is the eigenvector corresponding to the smallest eigenvalue of V T V. The eigenvector is normalized to obtain the required vector b, and then the matrix B is obtained. To obtain the matrix B, under the condition of a scale factor λ, the internal parameter is solved:

[0062]

[0063] where B 11 ,B 12 ,B 13 ,B 21 ,B 22 ,B 23 ,B 31 ,B 32 ,B 33 respectively represent the corresponding position elements in the matrix B, f x , f y , f s ,c x and c y are the elements in the internal parameter matrix K, f x , f y , f s represents a scaling ratio relationship, (c x ,c y), which represents the pixel coordinates of the camera principal point. After solving the intrinsic matrix K, the extrinsic matrix of each image is solved using the homography matrix H and the intrinsic matrix K of each image:

[0064]

[0065] In the formula, λ is the scale factor, r1, r2, and r3 are the three column vectors of the rotation matrix R, h1, h2, and h3 are the three column vectors of the homography matrix H, t represents the translation matrix M, and K is the intrinsic matrix;

[0066] Solve for the distortion parameters: Let (x s1 , y s1 ) be the ideal image point coordinates without considering distortion, and (x s , y s ) be the actual image point coordinates. Then, x ’ and y ’ are the abscissa and ordinate of the camera coordinate point after z - coordinate normalization, x ” and y ” are the abscissa and ordinate of the camera coordinate point after distortion, f x , f y , f s , c x and c y are the elements in the intrinsic parameter matrix K:

[0067]

[0068]

[0069] Subtracting gives:

[0070]

[0071] Arrange it into the matrix operation form regarding distortion:

[0072]

[0073]

[0074]

[0075] Among them, the X matrix is the distortion parameter, k1, k2, and k3 represent the radial distortion coefficients, p1 and p2 represent the tangential distortion coefficients; the D matrix is known; in the E matrix, (x s1 , y s1 ) are the ideal image point coordinates, (x s, y s ), which are the actual image point coordinates. The "T" in the upper right corner of the matrix indicates the transpose of the matrix. Therefore, the E matrix is also known;

[0076] When there are multiple points, use the matrix least squares algorithm to solve for the values of the camera distortion parameters k1, k2, k3, p1, and p2;

[0077] (31.4) Correct the original image using all the coefficients. The two corrected images lie in a plane and are parallel to each other;

[0078] (32) Generate a dense point cloud based on stereo matching and triangulation:

[0079] By performing feature matching and disparity calculation on the images captured by the left and right cameras, obtain the disparity value of each pixel point in the scene;

[0080] Utilize the disparity information and combine it with the parameters of the camera to calculate the coordinates of the object in three-dimensional space through triangulation;

[0081] (32.1) The left and right cameras simultaneously observe a point P in space w with coordinates (x w , y w , z w ). O l and O r are the optical centers of the left and right cameras respectively. P l and P r are the image point pairs of the spatial point P w after epipolar alignment in the left and right cameras. The coordinates of P l are (x l , y l ), and the coordinates of P r are (x r , y r ). Here, epipolar alignment is a method used in binocular stereo vision to increase the search efficiency of matching points;

[0082] Let the distance between the optical centers of the two cameras be b, the focal lengths of the cameras be f, and z be the depth information from the spatial point to the imaging plane. Then the relationship between P w and P l , P r is obtained according to the similarity relationship of triangles:

[0083]

[0084] Let d = x l - x r to get:

[0085]

[0086] When the disparity d of a certain point in the two cameras of the binocular stereo vision system is obtained, the depth information of this point can be obtained, that is, the three-dimensional coordinates are restored;

[0087] (33) Output the correspondence between the dense point cloud and the pixels of the left camera image;

[0088] Taking the left camera image as a template, the coordinates of the three-dimensional space points all have corresponding pixel coordinates in the image.

[0089] Further, based on the received left camera image, the structure from motion recovery algorithm is used to generate a sparse point cloud by multi-view triangulation, and the correspondence between the sparse point cloud and the pixels of the left camera image is output, as follows:

[0090] (41) Detect and match feature points and perform matching and geometric verification. Use the scale-invariant feature transform algorithm to extract features, represent the feature points by a 128-dimensional numerical descriptor, match the feature points, screen the matching point pairs based on the similarity of the descriptor vectors, and eliminate the abnormal matching point pairs through geometric verification;

[0091] (42) Use multi-view triangulation to generate a sparse point cloud. Initialize and select two matching images, set the pose of one image as the identity matrix, estimate the transformation matrix using the matching point pairs, decompose to obtain the pose of the other image, based on the camera poses of the two images, calculate the coordinates of the matching feature points in the three-dimensional space through triangulation to generate an initial sparse point cloud. As more images are added, match the new added image with the feature points in the existing point cloud, and iteratively optimize the poses of all images and the three-dimensional point coordinates to gradually improve the sparse reconstruction of the scene;

[0092] X W represents a three-dimensional point in space, O1 and O2 are the optical centers of the two cameras, and the plane formed by these two camera optical centers and point X W is called the epipolar plane, and I1 and I2 are the physical image planes corresponding to these two cameras;

[0093] In these two cameras, the pixel coordinates of point X W are x1 and x 2， respectively, and they represent two projection observations of the three-dimensional point X W ; x1', x2' are the normalized plane coordinates corresponding to x1, x2 respectively; the internal parameter matrices of the two cameras are K1, K2; the coordinate transformation matrix of the camera coordinate system is R t, and there are constraints:

[0094]

[0095]

[0096] Among them, E is the essential matrix and F is the fundamental matrix. After feature extraction and matching, the symbols in the upper right corner of the matrix represent matrix operations. According to the above constraints, the essential matrix and the fundamental matrix are calculated, and then the external parameters of the camera are obtained;

[0097] Let the point X W have pixel coordinates x1 = (u1, v1, 1) and x2 = (u2, v2, 1) as the feature points matched in two images, where u1, v1, u2, and v2 are the coordinate values corresponding to the x-axis and y-axis respectively, and we get:

[0098]

[0099] where F 11 , F 12 , F 13 , F 21 , F 22 , F 23 , F 31 , F 32 , F 33 respectively represent the corresponding position elements in the fundamental matrix F. Let f = (F 11 , F 12 , F 13 , F 21 , F 22 , F 23 , F 31 , F 32 , F 33 ) and we get:

[0100]

[0101] Each pair of matching points can provide a constraint for calculating the F matrix. Therefore, at least 8 pairs of matching feature points are required to solve for F. When the number of matching feature points is greater than 8, the least squares method or the RANSAC algorithm is used to solve the F matrix; Given the camera parameters and the matching points, the coordinates of the 3D points are obtained through triangulation. Let the projection matrix of the i-th camera be: P i = K i [R i , t i , P i is the projection matrix of the i-th camera, K i is the internal parameter matrix of the i-th camera, R i is the rotation matrix of the i-th camera, t i is the translation matrix of the i-th camera, the coordinates of the 3D point in space are X W = [x, y, z, 1] T , and the image coordinates of the corresponding feature point are x i = [u i , v i , 1]T , u i , v i is the abscissa and ordinate of point x i , and d i represents the depth of X W under i viewpoints, that is, the coordinate z in the camera coordinate system c Then:

[0102]

[0103] Cross multiply both sides by x i We get:

[0104]

[0105] Furthermore, we get:

[0106]

[0107] P i1 , P i2 , P i3 are the horizontal vectors at the corresponding positions of P i respectively. The rank of the system of equations is 2, and at least one pair of matching feature points is required to solve the coordinates of the three-dimensional point X;

[0108] (43) The correspondence between the sparse point cloud and the camera image pixels uses the input left camera image as a template, and the coordinates of the three-dimensional space points all have corresponding pixel coordinates in multiple images.

[0109] Furthermore, according to the correspondence between the obtained dense point cloud and the sparse point cloud and the image pixels respectively, calculate the correspondence between the dense point cloud and the sparse point cloud to obtain the preliminary alignment of the point cloud positions and complete the rough stitching of the point cloud, as follows:

[0110] (51) Each three-dimensional point in the dense point cloud can be in one-to-one correspondence with the pixel coordinates in the left camera image, and each three-dimensional point in the sparse point cloud can be in one-to-one correspondence with the pixel coordinates in the input image. By using the two-dimensional coordinates of the image pixels as the index value, the three-dimensional points in the dense point cloud and the three-dimensional points in the sparse point cloud can be associated, so as to obtain the spatial position correspondence between the dense point cloud and the sparse point cloud;

[0111] For each three-dimensional point X in the sparse point cloud xs , its coordinates are (x xs , y xs , z xs ), and there is a pixel point X t in the left camera image corresponding to it. The coordinates of X t are (x t , y t ), and the same pixel point X tIt must also correspond to a certain three-dimensional point X in the dense point cloud cm , X cm has coordinates (x cm , y cm , z cm )), then taking X t as the index, matching the corresponding three-dimensional points in the sparse point cloud and the dense point cloud, calculating the difference between the point cloud coordinates, and averaging the coordinate differences of the matching points to obtain the translation matrix M q , x q , y q , z q are the translation amounts corresponding to the x, y, and z directions respectively, where n represents the total number of pixel points X t , and t is the index number of the pixel point:

[0112] ;

[0113] where x xst , x cmt , y xst , y cmt , z xst , z cmt represent the t-th elements of x xs , x cm , y xs , y cm , z xs , z cm respectively;

[0114] (52) Based on the spatial position correspondence relationship between the dense point cloud and the sparse point cloud, calculate a global translation matrix to align the dense point cloud as a whole to the position of the sparse point cloud, and complete the rough stitching of the point cloud.

[0115] Furthermore, on the basis of the rough stitching of the point cloud, use the iterative closest point algorithm to perform precise stitching of the point cloud, minimizing the error in the overlapping area of the point cloud, as follows:

[0116] (61) Calculate the overlapping area of the point cloud, perform uniform sampling and match the closest points. Use the voxel hull represented by the low resolution of the volume occupied by the point cloud to determine the overlapping area. The object space is subdivided into a voxel structure, and points are uniformly selected in the object space so that the points in the overlapping area are evenly distributed. Match again, find the closest points of the selected subset of points to be matched in another point cloud, and reject incorrect correspondences based on the compatibility of the points;

[0117] Voxelize the point clouds P1 and P2, define the voxel resolution d v , and divide the object space V into a three-dimensional grid:

[0118]

[0119] Among them, v i，j，k represents the voxel located at the index (i, j, k), and the voxel center coordinate C i，j，k is:

[0120]

[0121] Point X v , with coordinates (x v , y v , z v ), and its corresponding voxel index is calculated as follows:

[0122]

[0123] Traverse all voxels. If there are points from both P1 and P2 in a voxel, then this voxel belongs to the overlapping region V overlap :

[0124]

[0125] Perform uniform sampling in V overlap to obtain a sparse but uniformly distributed point set S. Let P be the point set of the fixed point cloud. For each voxel v i，j，k ∈V overlap , uniformly select a representative point X s within this voxel to form the point set S:

[0126]

[0127] That is, select the centroid of all points within this voxel as the sampling point. For the point set S after uniform sampling, find the nearest point match in the point set P ’ of the other point cloud. For each sampling point X s , find the nearest neighbor point X ’ in P ’ s :

[0128]

[0129] Use the geometric compatibility of points to filter out incorrect matches. Let d s be the distance between the matching point pair (X s , X ’ S ), and set a threshold τ d for rejection:

[0130]

[0131] If d s > τ d, it is considered that the matching is incorrect and this point pair is removed;

[0132] (62) Based on the principle of minimizing the point-plane distance, iteratively calculate the optimal rotation and translation matrices until the convergence condition is met to complete the fine stitching:

[0133] Take the sparse point cloud as the target point cloud Q = {Y j = (x j , y j , z j ), T}, the dense point cloud as the point cloud to be registered P = {X i = (x i , y i , z i ), T}, the initial transformation parameters rotation matrix R o and translation matrix t o , the goal is to find R icp and t icp such that the transformed point cloud P ’ optimally matches Q:

[0134]

[0135] Estimate the normal vector N j of the target point cloud Q:

[0136]

[0137] where n xj , n xj , n xj are the corresponding coordinate values of the vector N j respectively, and the superscript T in the matrix represents the transpose operation of the matrix. For each point X i ∈ P, find the nearest point Y j in the target point cloud Q and obtain the normal vector N j . For each matching point pair (X i , Y j ), calculate the point-plane distance error d i :

[0138]

[0139] where N j T represents the transpose matrix of N j . Perform parameter optimization by minimizing the weighted sum of the squares of the point-plane distances of all corresponding points:

[0140]

[0141] Weight ωi Calculated from the median absolute deviation (MAD) of the residuals of the point cloud pairs to which the corresponding point pairs belong:

[0142]

[0143] where is the median, and σ k is the noise estimate of the k-th point cloud pair;

[0144] Let the parameter vector be β, and linearize the non-linear equation d i (x) = 0 + v i , where d i (x) is the point-plane distance error equation with respect to x, and v i represents the adjusted remaining distance, into a function related to x0: f Function:

[0145]

[0146] where d i (x) is the point-plane distance error equation with respect to x, vi represents the adjusted remaining distance; x0 represents the vector of approximate parameter values, A is the design matrix, which contains the partial derivatives of the equation with respect to the parameters at the point x0, v is the remaining distance, and △β is the parameter correction;

[0147] Solve for the parameter correction △β by iteration:

[0148]

[0149] In the formula, β0 is the initial value of the parameter, A is the design matrix, and the weight matrix P = diag(ω1, ω2…ω n ), and ω i is the corresponding weight;

[0150] The iteration update formula is:

[0151]

[0152] Converge when |△β| is below the set threshold;

[0153] The iterative optimization process includes external iteration and internal iteration. The external iteration includes re-establishing point cloud matching, calculating new normal vectors and point-plane distances, and filtering out incorrect matching points. The internal iteration includes linearized least squares solving for △β, updating the rotation parameters and translation parameters until convergence;

[0154] The point-to-plane distance between two corresponding points is defined as the orthogonal distance from a point to the fitted plane of another point. For the point clouds to be stitched, the transformation parameters are estimated by least squares adjustment, minimizing the sum of the squares of the point-to-plane distances, and finally, a transformation is performed, and the point clouds are transformed using the estimated parameters for point cloud stitching.

[0155] The present invention has at least the following beneficial effects:

[0156] 1. When the point cloud stitching of multi-view reconstruction is carried out in the present invention, the surface features of the underwater structure can be better restored without feature extraction, and it has strong accuracy and robustness, and high-precision point cloud stitching can be realized in the complex underwater environment. Therefore, it has good application prospects in the actual underwater structure detection and image stitching tasks.

[0157] 2. The present invention does not need to rely on feature markers, and the whole view reconstruction of the complex underwater structure can be completed in a single acquisition only through a binocular vision device.

[0158] 3. The present invention can better restore the three-dimensional topography of the concrete surface, and uses a stereo vision system to detect the underwater structure, which can improve the detection efficiency compared with the traditional manual underwater exploration.

[0159] Of course, it is not necessary for any product implementing the present invention to achieve all the above-mentioned advantages simultaneously. BRIEF DESCRIPTION OF THE DRAWINGS

[0160] Figure 1 is a schematic flow chart of the method of the present invention;

[0161] Figure 2 is a flow chart of dense point cloud reconstruction in the present invention;

[0162] Figure 3 is a schematic diagram of the binocular camera vision principle in the present invention;

[0163] Figure 4 is a flow chart of three-dimensional reconstruction of the structure from motion recovery algorithm in the present invention;

[0164] Figure 5 is a flow chart of rough stitching based on pixel index value in the present invention;

[0165] Figure 6 is a result diagram of fine stitching of multi-view underwater three-dimensional point clouds in the present invention;

[0166] Figure 7 is a schematic diagram of the result of fast stitching of the 3D reconstruction point cloud of the underwater structure surface in the present invention;

[0167] Figure 8 is a schematic diagram of the structure of the system of the present invention;

[0168] Figure 9 is a schematic diagram of the implementation process in the embodiment of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0169] Next, the technical solutions in the embodiments of the present disclosure will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present disclosure. Obviously, the described embodiments are only a part of the embodiments of the present disclosure, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments in the present disclosure without creative efforts shall fall within the protection scope of the present disclosure.

[0170] Embodiment 1:

[0171] The SIFT algorithm is a method that extracts stable, rotation and scale-invariant feature points from images through steps such as scale-space extreme value detection, key point localization, direction assignment, and descriptor generation, and is used for tasks such as image matching and object recognition.

[0172] Please refer to Figure 1 , the present invention provides a technical solution: a multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds, including:

[0173] S1. Receive the surface images of an underwater structure binocular camera, where the binocular camera includes a left camera and a right camera, generate a dense point cloud based on stereo matching and triangulation, and output the correspondence between the dense point cloud and the pixels of the left camera image, specifically as follows:

[0174] (S11) Use a stereo vision system with a resolution of 5472×3648 pixels to collect the surface images of the underwater structure. The system is calibrated using the standard checkerboard calibration method to obtain the internal and external parameters of the camera;

[0175] (S11.1) Camera calibration is a process of establishing the mutual position correspondence between the pixel coordinate system and the object points in space, and is also a process of determining the internal parameter matrix K, external parameter matrix Rt, and distortion coefficients of the camera;

[0176] The internal parameter matrix K of the camera contains the basic attributes of the camera and is used to describe the optical characteristics of the camera. Its form is:

[0177]

[0178] Where: f x , f y represents a scaling ratio relationship, which is calculated through the image distance f , pixel size s and the true angle θ of the pixel arrangement on the camera plane:

[0179]

[0180] Where: s x , sy is the width and height of the pixel size; (c x , c y ) represents the pixel coordinates of the camera principal point;

[0181] f s It is calculated from the image distance during imaging f , the horizontal pixel size s x and the true angle θ of the pixel arrangement in the sensor plane:

[0182]

[0183] △r is the mathematical model of radial distortion and is expressed as:

[0184]

[0185] The radial distortion is decomposed into the x and y axis directions to obtain the Gaussian radial distortion model. Among them, r is the distance from the image point to the center of the image plane, r 2 = x 2 + y 2 , k1, k2, and k3 represent the radial distortion coefficients, x and y respectively represent the abscissa and ordinate of the point in the ideal image coordinate system, and △X r represents the Gaussian radial distortion decomposed into the x axis direction, and △Y r in the x axis direction are respectively expressed as:

[0186]

[0187] P(r) is the mathematical model of tangential distortion. Among them, r is the distance from the image point to the center of the image plane, r 2 = x 2 + y 2 , p1 and p2 represent the tangential distortion coefficients and are expressed as:

[0188]

[0189] The tangential distortion is decomposed into the x and y axis directions. x and y respectively represent the abscissa and ordinate of the point in the ideal image coordinate system. △x d represents the tangential distortion decomposed into the x axis direction, and △y d represents the tangential distortion decomposed into the y axis direction and is respectively expressed as:

[0190]

[0191] In addition, the external camera parameters include a 3*3 rotation matrix R and a 3*1 translation matrix M to complete the world coordinate system ( O w -x w y w z w ) Transformation to the camera coordinate system ( O -[[]] xyz ) Affine transformation, where x w , y w , z w respectively represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the point ( x w , y w , z w ) in the three-dimensional world coordinate system, x , y , z respectively represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the point ( x , y , z ) in the three-dimensional camera coordinate system. The affine transformation is expressed as follows:

[0192]

[0193]

[0194] The matrix A is composed of the rotation matrix R and the translation matrix M, where R 11 , R 12 , R 13 , R 21 , R 22 , R 23 , R 31 , R 32 , R 33 respectively represent the corresponding position elements in the rotation matrix R, and M x , M y , M z respectively represent the corresponding position elements in the translation matrix M;

[0195] And in the second step, the camera coordinate system needs to be transformed to the image coordinate system ( O i – x i y i ), where x i , y i respectively represent the point ( xi , y i ), the x-axis coordinate and y-axis coordinate, and the coordinate transformation is expressed as:

[0196]

[0197] In the formula, f is the image distance during imaging, α is the scale factor of the coordinate, and the B matrix is composed of the image distance f Finally, the image coordinates are converted to the pixel coordinate system ( O s – x s y s ), where x s , y s respectively represent the x-axis coordinate and y-axis coordinate of the point ( x s , y s ) in the pixel coordinate system:

[0198]

[0199] In the formula, s x , s y are the width and height of the pixel size; (c x , c y ) represents the pixel coordinates of the camera principal point, θ is the true angle of the pixel arrangement on the camera plane, and the C matrix is composed of 1 / s x , -cotθ / s y , c x , 1 / ( s y sinθ), c y Combining and substituting all coordinate systems to obtain the camera imaging mathematical model:

[0200]

[0201]

[0202] (S11.2) The internal parameter matrix K and distortion coefficients of the left camera and the right camera correct the original images, and the two corrected images are located in a plane and parallel to each other;

[0203] For the camera coordinate system point (I, J, K), there is:

[0204]

[0205]

[0206]

[0207] Among them, k1, k2, and k3 represent radial distortion coefficients, and p1, p2 represent tangential distortion coefficients, where R 11 , R 12 , R 13 , R 21 , R 22 , R 23 , R 31 , R 32 , R 33 respectively represent the corresponding position elements in the rotation matrix R, M x , M y , M z respectively represent the corresponding position elements in the translation matrix M, r is the distance from the image point to the center of the image plane, r 2 = x 2 + y 2 , f x , f y , f s , c x and c y are the elements in the intrinsic parameter matrix K, x ’ and y ’ are the abscissa and ordinate of the camera coordinate point after z - coordinate normalization, x ” and y ” are the abscissa and ordinate of the distorted camera coordinate point, (x s , y s ) is the pixel coordinate of the image point. The intrinsic and extrinsic parameters are calibrated using a checkerboard calibration board. For each feature point on the calibration board, the following equation is satisfied:

[0208]

[0209] In the above equation, r1, r2, r3 are the three column vectors of the rotation matrix R, t represents the translation matrix M, K is the intrinsic parameter matrix, α is the scale factor of the coordinate, (x w , y w , z w ) is the coordinate of the point on the checkerboard in the world coordinate system. Establish the homography mapping relationship between the calibration board point (x w , y w , z w ) and the image point (x s , y s ):

[0210]

[0211] Among them, H is the homography matrix, h 11 , h 12 , h 13 , h 21 , h 22 , h 23 , h 31 , h 32 respectively represent the corresponding position elements in the homography matrix H, s is an unknown scale factor, and H ’ =[h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 T The matrix operation equation of:

[0212]

[0213] As can be seen from the above, the equations of multiple corresponding points can be superimposed and written in the form of a non - homogeneous system of equations: AX = C; then use the least - squares method to solve this system of equations: X=(A T A) -1 A T C. At least 4 pairs of known spatial points and image points in a single photo are required to solve the homography matrix H. 8 equations are used to solve 8 unknowns, and each photo has a homography matrix;

[0214] And since , where h1, h2, h3 are the three column vectors of the homography matrix H, λ is the scale factor, r1, r2 are the column vectors of the first two columns of the rotation matrix R, and t is the translation matrix T. Since the rotation matrix is an orthogonal matrix, its rows or column vectors are pairwise orthogonal and unit vectors, thus obtaining two constraint equations:

[0215]

[0216] K is the camera intrinsic matrix, h1 T , h2 T respectively represent the transposes of the column vectors h1 and h2, K -T represents the transpose of the inverse matrix of the K vector, K -1 represents the inverse matrix of the K vector. Let:

[0217]

[0218] B is a symmetric matrix, B 11 , B​12 , B 13 , B 21 , B 22 , B 23 , B 31 , B 32 , B 33 Represent the corresponding position elements in matrix B respectively, and use a 6-dimensional vector to define the b matrix. The T in the upper right corner of the matrix represents the transpose of the matrix:

[0219]

[0220] It can be seen that:

[0221]

[0222] h1 and h2 are the first two column vectors of the homography matrix H, h i You can choose h1, h2, h according to the value of i j You can choose h1, h2, h according to the value of j i T h i Inverse matrix, where v ij T The vector is represented as follows, and the T in the upper right corner of the matrix represents the transpose of the matrix:

[0223]

[0224] h 11 , h 12 , h 13 , h 21 , h 22 , h 23 Represent the corresponding position elements in the homography matrix H, then the constraint equation is expressed by the matrix v ij and b are written in matrix operation form:

[0225]

[0226] Among them, v 11 T Indicates v ij T Corresponding position element, v 12 , v 22 , respectively representing v ij Corresponding position elements, from the above formula we can see that multiple groups of v ij Substitute the expression to write it in the form of a homogeneous equation solution: Vb=0; the least squares solution of the homogeneous equation system is V TThe eigenvector corresponding to the minimum eigenvalue of V. Normalize the eigenvector to obtain the required vector b, and then obtain matrix B. There are 6 unknowns in vector b, and at least 3 sets of v are required. ij That is, 3 H matrices and 6 equations to solve 6 unknowns. Therefore, at least 3 different calibration pictures are required to obtain matrix b. After obtaining matrix b, construct matrix B. Under the condition of a scale factor difference (B = λK -T K -1 , where λ is the scale factor), solve the internal parameter:

[0227]

[0228] where B 11 , B 12 , B 13 , B 21 , B 22 , B 23 , B 31 , B 32 , B 33 respectively represent the corresponding elements in matrix B. f x , f y , f s , c x and c y are the elements in the internal parameter matrix K. f x , f y , f s represents a scaling ratio relationship. (c x , c y ) represents the pixel coordinates of the camera principal point. After solving the internal parameter matrix K, the external parameter matrix of each image can be solved using the homography matrix H and the internal parameter matrix K of each image:

[0229]

[0230] where λ is the scale factor, r1, r2, r3 are the three column vectors of the rotation matrix R, h1, h2, h3 are the three column vectors of the homography matrix H, t represents the translation matrix M, and K is the internal parameter matrix. Next, solve the distortion parameters: Let (x s1 , y s1 ) be the ideal image point coordinates without considering distortion, (x s , y s ) be the actual image point coordinates, x ’ and y ’ are the abscissa and ordinate of the camera coordinate point after z coordinate normalization, x” and y ” are the abscissa and ordinate of the camera coordinate points after distortion, f x , f y , f s , c x and c y are the elements in the internal parameter matrix K, then:

[0231]

[0232]

[0233] Subtracting gives:

[0234]

[0235] Rearranging into a matrix operation form regarding distortion:

[0236]

[0237]

[0238]

[0239] Among them, the X matrix is the distortion parameter to be solved, k1, k2, k3 represent the radial distortion coefficients, and p1, p2 represent the tangential distortion coefficients; the camera internal and external parameters required in the D matrix have been obtained before, so the D matrix is known; in the E matrix, (x s1 , y s1 ) are the ideal image point coordinates, (x s , y s ) are the actual image point coordinates, and the T in the upper right corner of the matrix represents the transpose of the matrix, so the E matrix is also known; when there are multiple points, it can be rearranged into a matrix operation form, and the values of the camera distortion parameters k1, k2, k3, p1, p2 can be solved using the matrix least squares algorithm;

[0240] Finally, use all the coefficients to correct the original image, and the two corrected images are in a plane and parallel to each other;

[0241] (S12) Under the condition of ensuring image field of view overlap, the system translates along the surface of the structure to collect multiple groups of images, generates a dense point cloud based on stereo matching and triangulation, and outputs the correspondence between the dense point cloud and the pixels of the left camera image, specifically as follows:

[0242] By performing feature matching and disparity calculation on the images captured by the left and right cameras, the disparity value of each pixel point in the scene is obtained;

[0243] Using the disparity information and combining it with the parameters of the camera, the coordinates of the object in the three-dimensional space are calculated through triangulation;

[0244] The left and right cameras simultaneously observe a point P w (x w , y w , z w ), O l and O r are the optical centers of the left and right cameras respectively. P l (x l , y l ) and P r (x r , y r ) are the image point pairs of the spatial point P w after epipolar alignment in the left and right cameras respectively. Among them, epipolar alignment is a method used in binocular stereo vision to increase the search efficiency of matching points. Assuming that the distance between the optical centers of the two cameras is b and the focal lengths of the cameras are both f, then P w and P l , P r 's relationship can be obtained according to the similarity relationship of triangles:

[0245]

[0246] Let d = x l -x r It can be obtained:

[0247]

[0248] After obtaining the disparity d of a certain point in the binocular stereo vision system's two cameras, the depth information of this point can be obtained, that is, the three-dimensional coordinates are restored;

[0249] (S13) Output the correspondence between the dense point cloud and the pixels of the left camera image;

[0250] Taking the left camera image as a template, the coordinates of the three-dimensional space points all have corresponding pixel coordinates in the image;

[0251] Specifically, the specific acquisition device for the underwater stereo vision image is an underwater binocular camera. The camera optical method is used to obtain two-dimensional images with overlapping fields of view of the underwater structure and its surface at different angles. The binocular camera includes a left camera and a right camera, and the working distance between the two cameras is about 300 mm;

[0252] Regarding the technical solution of this embodiment, such as Figure 2The following is a flowchart of dense point cloud reconstruction. It mainly uses the binocular vision principle and two images with overlapping areas collected by a binocular camera to reconstruct the underwater structure surface topography. The main processes include underwater calibration of the calibration board, selection of the image overlapping area, and calculation of dense point cloud coordinates. For the binocular system calibration using the calibration board, the calibration board and the binocular system are both placed in the underwater environment. Calibration images of the calibration board are collected at different object distances, and the internal and external parameters of the left and right cameras are calculated. Based on the internal and external parameters, epipolar correction and distortion correction are performed. Then, the left and right camera images of the binocular camera are input. By performing feature matching and disparity calculation on the images taken by the left and right cameras, the disparity value of each pixel point in the scene is obtained. Using the disparity information and combining with the camera parameters, the coordinates of the object in the three-dimensional space are calculated through triangulation. Taking the left camera image as a template, the corresponding relationship between the dense point cloud and the pixels of the left camera image is output;

[0253] Regarding the technical solution of this embodiment, as Figure 3 shown in the binocular camera vision principle diagram, O c –x c y c z c and O I –x1y1z1 represent the camera coordinate system and pixel coordinate system of each camera respectively. P1(x1, y1, z1) and P2(x2, y2, z2) are the corresponding imaging points of point P on the imaging planes of camera A and camera B respectively. The optical center of camera A (O c1 ), imaging point P1, and actual point P are collinear. Similarly, the optical center of camera B (O c2 ), imaging point P2, and actual point P are also collinear. Therefore, by calculating the intersection coordinates of the lines O c1 P1 and O c2 P2, the coordinates of point P in the three-dimensional space can be obtained; the coordinates of all point clouds are based on the coordinate system of the left camera, and the two-dimensional points on the left camera image are converted into three-dimensional points in space through matrix transformation;

[0254] S2. Based on the received left camera image, the structure from motion recovery algorithm is used to generate a sparse point cloud using multi-view triangulation and output the corresponding relationship between the sparse point cloud and the pixels of the left camera image, specifically as follows:

[0255] (S21) Detect and match feature points and perform matching and geometric verification. The scale-invariant feature transform algorithm is used to extract features, and the feature points are represented by a 128-dimensional numerical descriptor. The feature points are matched, and the matching point pairs are screened based on the similarity of the descriptor vectors, and the abnormal matching point pairs are removed through geometric verification;

[0256] (S22) Generate a sparse point cloud using multi-view triangulation. Initialize by selecting two matching images, set the pose of one image as the identity matrix, estimate the transformation matrix using the matching point pairs, decompose to obtain the pose of the other image, and based on the camera poses of the two images, calculate the coordinates of the matching feature points in 3D space through triangulation to generate an initial sparse point cloud. As more images are added, match the new image with the feature points in the existing point cloud, and iteratively optimize the poses of all images and the 3D point coordinates to gradually improve the sparse reconstruction of the scene;

[0257] X W represents a 3D point in space, and O₁ and O₂ are the optical centers of two cameras; the plane formed by these two camera optical centers and point X W is called the epipolar plane; I₁ and I₂ are the physical image planes corresponding to these two cameras. In these two cameras, the pixel coordinates of point X W are x₁ and x 2， respectively, and they represent two projection observations of the 3D point X W ; x₁', x₂' are the normalized plane coordinates corresponding to x₁, x₂ respectively; the internal parameter matrices of the two cameras are k₁, k₂; the coordinate transformation matrix of the camera coordinate system is R t, and there are constraints:

[0258]

[0259]

[0260] where E is the essential matrix and F is the fundamental matrix; after feature extraction and matching, the symbols in the upper right corner of the matrix represent matrix operations, and the essential matrix and fundamental matrix can be calculated according to the above constraint conditions, and then the external parameter information of the camera can be obtained;

[0261] Let the pixel coordinates of point X W be x₁ = (u₁, v₁, 1), x₂ = (u₂, v₂, 1) as the matching feature points of two images, where u₁, v₁, u₂, v₂ are the coordinate values corresponding to the x-axis and y-axis respectively, and we get:

[0262]

[0263] where F 11 ,F 12 ,F 13 ,F 21 ,F 22 ,F 23 ,F 31 ,F 32 ,F 33 represent the corresponding position elements in the fundamental matrix F respectively. Let f = (F 11 ,F 12 ,F13 , F 21 , F 22 , F 23 , F 31 , F 32 , F 33 )have to:

[0264]

[0265] Each pair of matching points can provide a constraint for calculating the F matrix, so solving F requires at least 8 pairs of matching feature points. When the number of matching feature points is greater than 8, the least squares method or RANSAC algorithm is used to solve the F matrix. If the camera parameters and matching points are known, the coordinates of the three-dimensional points are obtained by triangulation. Let the projection matrix of the i-th camera be: P i =K i [R i ,t i ], P i is the projection matrix of the i-th camera, K i is the intrinsic parameter matrix of the i-th camera, R i is the rotation matrix of the i-th camera, t i is the translation matrix of the i-th camera, and the coordinates of the three-dimensional point in space are X W =[x,y,z,1] T , the image coordinates of the corresponding feature points are x i =[u i ,v i ,1] T ,u i , v i is x i The horizontal and vertical coordinates of the point, d i Represents X W The depth at the i-th perspective, that is, the coordinate z in the camera coordinate system c but:

[0266]

[0267] Cross product x on both sides i We can get:

[0268]

[0269] Further we can get:

[0270]

[0271] P i1 , P i2 , P i3 P iThe horizontal vectors at the corresponding positions. The rank of this system of equations is 2. At least one pair of matching feature points is required to solve for the coordinates of the 3D point X. In fact, for a single 3D point, there are multiple feature matching point pairs in the input images from multiple perspectives. The coordinates of the 3D point X can be accurately solved by the least squares method or the RANSAC method;

[0272] (S23)The correspondence between the sparse point cloud and the pixels of the camera image. Taking the input left camera image as a template, the coordinates of the 3D space points all have corresponding pixel coordinates in multiple images;

[0273] Specifically, as Figure 4 shown in the flowchart of 3D reconstruction of the Structure from Motion (SFM) algorithm: First, import the images of the left camera collected by the binocular camera system, which are image sequences taken at different angles and time points. During the process of retrieving matches, use the SIFT feature extraction method to find sparse feature points in the images, and use a numerical descriptor (128-dimensional vector) to describe the feature points, then match the image pairs again. Secondly, match the feature points on the matched images, use the similarity measure of the 128-dimensional descriptor vectors of the feature points to perform feature point matching, and geometrically verify the matched points. After passing the verification, initialize and select two matched images, set the pose of one of the images as the identity matrix, and then estimate the matrix through the matched point pairs between them. Decompose the matrix to obtain the pose of the next image (select and reconstruct the next image). After estimating the poses of the two images, calculate the positions of the matched feature points in the 3D space through triangulation, and then perform bundle adjustment and filter out outliers to generate an initial 3D sparse point cloud. With the addition of more images, continuously match the new images with the feature points in the existing point cloud to gradually optimize the poses of all images and the 3D point coordinates. Finally, export the correspondence between the point cloud and the image pixels simultaneously;

[0274] S3. According to the correspondences between the obtained dense point cloud and the sparse point cloud and the image pixels respectively, calculate the correspondence between the dense point cloud and the sparse point cloud to obtain a preliminary alignment of the point cloud positions and complete the rough stitching of the point cloud, as follows:

[0275] (S31)Each 3D point in the dense point cloud can be in one-to-one correspondence with the pixel coordinates in the left camera image, and each 3D point in the sparse point cloud can be in one-to-one correspondence with the pixel coordinates in the input image. By using the two-dimensional coordinates of the image pixels as index values, the 3D points in the dense point cloud and the 3D points in the sparse point cloud can be associated, thereby obtaining the spatial position correspondence between the dense point cloud and the sparse point cloud;

[0276] For each 3D point X xs (x xs ,y xs ,z xs ) in the sparse point cloud, there is a pixel point X in the left camera imaget (x t , y t ), there is a corresponding one, and similarly, this pixel point X t must also correspond to a certain three-dimensional point X cm (x cm , y cm , z cm ) in the dense point cloud. Then, using X t as an index, match the corresponding three-dimensional points in the sparse point cloud and the dense point cloud, calculate the difference between the point cloud coordinates, and take the average of the coordinate differences of the matching points to obtain the translation matrix M q , where n represents the total number of pixel points X t , and t is the index number of the pixel point:

[0277]

[0278] where x xst , x cmt , y xst , y cmt , z xst , z cmt respectively represent the t-th elements of x xs , x cm , y xs , y cm , z xs , z cm ;

[0279] (S32) Based on the spatial position correspondence between the dense point cloud and the sparse point cloud, calculate a global translation matrix to align the dense point cloud as a whole to the position of the sparse point cloud, and complete the rough stitching of the point cloud;

[0280] Specifically, as Figure 5 shown is the flowchart of rough stitching based on the pixel index value. In the initial state, there are phenomena of coordinate mismatch and excessive distance between different reconstructed dense point clouds and the SFM sparse point cloud. Using the correspondence between the sparse point cloud and image pixels and the correspondence between the dense point cloud and image pixels, taking the left camera image pixels as the index value, perform point cloud to pixel matching, calculate the translation value of each binocular point cloud, and perform rough stitching by translating the binocular point cloud to align the dense point cloud as a whole to the position of the sparse point cloud, ensuring the approximate coincidence between the point clouds;

[0281] S4. On the basis of the rough stitching of the point cloud, use the iterative closest point algorithm to perform precise stitching of the point cloud and minimize the error in the overlapping area of the point cloud, as follows:

[0282] (S41) Calculate the overlapping region of the point clouds, perform uniform sampling and match the nearest points. Determine the overlapping region using a voxel hull represented by the low resolution of the volume occupied by the point cloud. The object space is subdivided into a voxel structure, and points are uniformly selected in the object space such that the points within the overlapping region are evenly distributed. Then perform the matching again, find the nearest points of the selected subset of points to be matched in the other point cloud, and reject incorrect correspondences based on the compatibility of the points.

[0283] Voxelize the point clouds P1 and P2, and define the voxel resolution d v , divide the object space V into a three-dimensional grid:

[0284] [[ID=�10]]

[0285] where, v i，j，k represents the voxel located at the index (i, j, k), and the voxel center coordinate C i，j，k is:

[0286]

[0287] For the point X v (x v , y v , z v ), the calculation of the voxel index it belongs to is as follows:

[0288]

[0289] Traverse all voxels. If there are points from both P1 and P2 within a voxel, then this voxel belongs to the overlapping region V overlap [[ID=�37]]:

[0290] )

[0291] Perform uniform sampling in V overlap to obtain a sparse but evenly distributed point set S. Let P be the point set of the fixed point cloud. For each voxel v i，j，k ∈ V overlap , uniformly select a representative point X s within this voxel to form the point set S:

[0292]

[0293] That is, select the centroid of all points within this voxel as the sampling point. For the uniformly sampled point set S, find the nearest point match in the point set P ’ of the other point cloud. For each sampling point X s , find the nearest neighbor point X ’ in P ’ s :

[0294]

[0295] To improve the accuracy of matching, the geometric compatibility of points is used to filter out incorrect matches. Let d s be the distance between the matching point pair (X s , X ’ S ). Set a threshold τ d for rejection:

[0296]

[0297] If d s > τ d , the matching is considered incorrect and the point pair is rejected;

[0298] (S42) Based on the principle of minimizing the point-plane distance, iteratively calculate the optimal rotation and translation matrices until the convergence condition is met to complete the fine stitching:

[0299] Take the sparse point cloud as the target point cloud Q = {Y j = (x j , y j , z j ) T}, the dense point cloud as the point cloud to be registered P = {X i = (x i , y i , z i ) T}, the initial transformation parameters rotation matrix R o and translation matrix t o . The goal is to find R icp and t icp such that the transformed point cloud P ’ is the optimal match for Q:

[0300]

[0301] Estimate the normal vector N j of the target point cloud Q:

[0302]

[0303] It can be calculated by PCA or obtained from an existing model, where n xj , n xj , n xj are the corresponding coordinate values of the vector N j respectively. The symbol in the upper right corner of the matrix represents the transpose operation of the matrix. For each point X i ∈ P, find the nearest point Y j in the target point cloud Q and obtain the normal vector N j . For each matching point pair (Xi , Y j ), calculate the point-plane distance error d i :

[0304]

[0305] where N j T denotes the transpose matrix of N j , and parameter optimization is performed by minimizing the weighted sum of squares of point-plane distances for all corresponding points:

[0306]

[0307] The weight ω i is calculated from the median absolute deviation (MAD) of the residuals of the point cloud pair to which the corresponding point pair belongs:

[0308]

[0309] where is the median, and σ k is the noise estimate of the k-th point cloud pair;

[0310] Let the parameter vector be β, and linearize the non-linear equation d i (x) = 0 + v i , where d i (x) is the point-plane distance error equation with respect to x, and v i represents the adjusted remaining distance, into a function related to x0: f Function:

[0311]

[0312] where x0 represents the vector of approximate parameter values, A is the design matrix, which contains the partial derivatives of the equation with respect to the parameters at the point x0, v is the remaining distance, and △β is the parameter correction;

[0313] Solve for the parameter correction △β by iteration:

[0314]

[0315] In the formula, β0 is the initial value of the parameter, A is the design matrix, and the weight matrix P = diag(ω1, ω2…ω n ), and ω i is the corresponding weight;

[0316] The iteration update formula is:

[0317]

[0318] Until convergence, i.e., |△β| is lower than the set threshold; the iterative optimization process includes external iteration and internal iteration. The external iteration includes re - establishing point cloud matching, calculating new normal vectors and point - plane distances, and filtering out incorrect matching points. The internal iteration includes linearizing the least - squares solution for △β, updating the rotation parameters and translation parameters until convergence;

[0319] The point - to - plane distance between two corresponding points is defined as the orthogonal distance from a point to the fitted plane of another point. For the point clouds to be stitched, the transformation parameters are estimated through least - squares adjustment, minimizing the sum of the squares of the point - to - plane distances, and finally performing transformation, using the estimated parameters to transform the point clouds for point cloud stitching;

[0320] It can improve the orientation of a single loose point cloud relative to a single fixed point cloud, and it is required that the point clouds are already roughly aligned;

[0321] Specifically, as Figure 6 is the fine - stitching result diagram of multi - view underwater three - dimensional point clouds:

[0322] Initial state: The point clouds are already roughly aligned. Use the voxel hull represented by the low resolution of the volume occupied by the point clouds to find and determine the overlapping area (finding the overlapping area). The object space is subdivided into a voxel structure, and points are uniformly selected in the object space (selection of points in the first point cloud) so that the points in the overlapping area are evenly distributed. Then perform matching again, find the nearest points of the selected subset in another point cloud (matching of selected points and nearest points), and reject incorrect correspondences based on the compatibility of the points;

[0323] The point - to - plane distance between two corresponding points is defined as the orthogonal distance from a point to the fitted plane of another point. For the point clouds to be stitched, the transformation parameters are estimated through least - squares adjustment, which will minimize the sum of the squares of the point - to - plane distances, and finally perform transformation, using the estimated parameters to transform the point clouds (minimization of point - to - plane distance);

[0324] Use the estimated parameters to transform the point clouds for point cloud stitching, and finally test the convergence criterion. If not met, the process will start over until the result converges, and export the best rotation and translation matrices (fine - stitching result).

[0325] Next, the present invention will be further elaborated in combination with specific implementation cases:

[0326] Such as Figure 7This is an exemplary embodiment of the rapid stitching result of the underwater structure surface 3D reconstruction point cloud using the method proposed in the present invention, showing the underwater reconstruction effect at different camera angles. The four sub-figures respectively correspond to the reconstruction results at tilt angles of 25°, 20°, 15°, and 10°. Each figure presents the three-dimensional point cloud model after rapid stitching of the point cloud. The coordinate axes represent the X, Y, and Z directions respectively, and the Z-axis uses color mapping to display the depth distribution. It can be observed that as the tilt angle decreases, the spatial distribution and depth information of the reconstruction area gradually change, and the restoration degree and detail performance of the surface structure increase with the decrease of the angle, reflecting the influence of the camera pose on the underwater structure three-dimensional point cloud stitching and reconstruction effect. The method proposed in this embodiment can better restore the surface features of the underwater structure when stitching the point clouds of multi-view reconstruction, without the need for feature extraction, and has strong accuracy and robustness, showing good application prospects in actual underwater structure detection and image stitching tasks.

[0327] Figure 9 As shown in the schematic diagram of the implementation process in the embodiment, the test environment is constructed inside a transparent water tank, and the underwater experimental conditions are simulated by injecting water. A concrete sample with a specification of 40 cm (length) × 80 cm (width) is set at the center of the water tank as the observation object. A binocular vision acquisition system is used for detection, which includes two symmetrically distributed industrial cameras (marked as Camera A and Camera B), both of which are fixed inside the water tank to collect data on the surface of the sample. During the test process, the system acquires data according to the preset three-dimensional acquisition path, and the entire system is connected to an external control terminal through a waterproof cable to achieve image acquisition and data transmission in the underwater environment.

[0328] Embodiment Two:

[0329] As Figure 8 shown, according to the second aspect of the present invention, the present invention provides a multi-view underwater three-dimensional point cloud stitching system based on the fusion of sparse and dense point clouds, which is used to implement the above-mentioned multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds, and includes:

[0330] A dense point cloud reconstruction module, which is used to receive the surface images of the underwater structure binocular cameras. The binocular cameras include a left camera and a right camera, generate a dense point cloud based on stereo matching and triangulation, and output the corresponding relationship between the dense point cloud and the pixels of the left camera image;

[0331] A sparse point cloud reconstruction module, which is used to generate a sparse point cloud based on the received left camera image by using the structure from motion algorithm and multi-view triangulation, and output the corresponding relationship between the sparse point cloud and the pixels of the left camera image;

[0332] A rough stitching module, which is used to calculate the correspondence between the dense point cloud and the sparse point cloud according to the correspondence between the obtained dense point cloud and sparse point cloud and the image pixels, obtain the preliminary alignment of the point cloud positions, and complete the rough stitching of the point cloud;

[0333] A fine stitching module, which, based on the rough stitching of the point cloud, uses the iterative closest point algorithm to accurately stitch the point cloud and minimize the error in the overlapping area of the point cloud.

[0334] Embodiment III:

[0335] The present invention provides a terminal device, including a memory, a processor, and a computer program stored in the memory and capable of running on the processor. The memory stores a computer program capable of running on the processor. When the processor loads and executes the computer program, the above-mentioned multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds is adopted.

[0336] Embodiment IV:

[0337] The present invention provides a storage medium containing computer-executable instructions, and the computer-executable instructions are used to execute the above-mentioned multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds when executed by a computer processor.

Claims

1. A multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds, characterized in that It includes the following steps: Receive the surface images of the binocular cameras of the underwater structure. The binocular cameras include a left camera and a right camera. Generate a dense point cloud based on stereo matching and triangulation, and output the correspondence between the dense point cloud and the pixels of the left camera image. Based on the received left camera image, use the structure from motion algorithm to generate a sparse point cloud by multi-view triangulation, and output the correspondence between the sparse point cloud and the pixels of the left camera image. According to the correspondences between the obtained dense point cloud and sparse point cloud and the image pixels respectively, calculate the correspondence between the dense point cloud and the sparse point cloud to obtain a preliminary alignment of the point cloud positions and complete the rough stitching of the point cloud. On the basis of the rough stitching of the point cloud, use the iterative closest point algorithm to perform precise stitching of the point cloud to minimize the error in the overlapping area of the point cloud. On the basis of the rough stitching of the point cloud, use the iterative closest point algorithm to perform precise stitching of the point cloud to minimize the error in the overlapping area of the point cloud, specifically as follows: Calculate the overlapping area of the point cloud, perform uniform sampling and match the closest points. Use the voxel shell represented by the low resolution of the volume occupied by the point cloud to determine the overlapping area. The object space is subdivided into a voxel structure, and points are uniformly selected in the object space so that the points in the overlapping area are evenly distributed. Match again, find the closest points of the selected subset of points to be matched in another point cloud, and reject incorrect correspondences according to the compatibility of the points. Voxelize the point clouds P1 and P2, and define the voxel resolution d v , and divide the object space V into a three-dimensional grid: where v i，j，k represents the voxel located at the index (i, j, k), and the voxel center coordinate C i，j，k is: Point X v , with coordinates (x v , y v , z v ), and its voxel index calculation is as follows: Traverse all voxels. If there are points from both P1 and P2 in a voxel, then this voxel belongs to the overlapping region V overlap : Uniform sampling is performed in V overlap to obtain a sparse but uniformly distributed point set S. Let P be the point set of the fixed point cloud. For each voxel v i，j，k ∈V overlap , a representative point X is uniformly selected within this voxel s to form the point set S: That is, the centroid of all points within this voxel is selected as the sampling point. For the point set S after uniform sampling, the nearest point matching is searched in the point set P of another point cloud. ’ For each sampling point X s , the nearest neighbor point X ’ is searched in P ’ s : Use the geometric compatibility of points to filter out incorrect matches. Let d s be the distance of the matching point pair (X s , X ’ S ). Set a threshold τ d for rejection: If d s > τ d , it is considered that the matching is incorrect and this point pair is excluded; Based on the principle of minimizing the point-to-plane distance, iteratively calculate the optimal rotation and translation matrices until the convergence condition is met to complete the fine stitching: Take the sparse point cloud as the target point cloud Q = {Y j = (x j , y j , z j ) T}, the dense point cloud as the point cloud to be registered P = {X i = (x i , y i , z i ) T}, the initial transformation parameters rotation matrix R o and translation matrix t o , the goal is to find R icp and t icp such that the transformed point cloud P ’ optimally matches Q: Estimate the normal vector N of the target point cloud Q j : where n xj ,n xj ,n xj are the corresponding coordinate values of the vector N j respectively. The superscript T in the upper right corner of the matrix represents the transpose operation of the matrix. For each point X i ∈P, find the nearest point Y j in the target point cloud Q, and obtain the normal vector N j . For each pair of matching points (X i , Y j ), calculate the point-plane distance error d i : where N j T denotes the transpose matrix of N j , and parameter optimization is performed by minimizing the sum of the weighted point-plane distance squares of all corresponding points: Weight ω i Calculated from the median absolute deviation (MAD) of the residuals of the point cloud pairs to which the corresponding point pairs belong: where is the median, and σ k is the noise estimate of the k-th point cloud pair; Let the parameter vector be β, and linearize the non - linear equation d i (x) = 0 + v i into a function related to x0 f : where d i (x) is the point-plane distance error equation with respect to x, vi represents the adjusted remaining distance; x0 represents the vector of approximate parameter values, A is the design matrix, containing the partial derivatives of the equation with respect to the parameters at the point x0, v is the remaining distance, and △β is the parameter correction amount; Iteratively solve the parameter correction amount △β: where β0 is the initial value of the parameter, A is the design matrix, and the weight matrix P = diag(ω1, ω2…ω n ), ω n is the corresponding weight; The iterative update formula is: Converge when |△β| is lower than the set threshold; The iterative optimization process includes external iteration and internal iteration. The external iteration includes re-establishing point cloud matching, calculating new normal vectors and point-to-plane distances, and filtering incorrect matching points. The internal iteration includes linearly solving the least squares for △β and updating the rotation parameters and translation parameters until convergence. The point-to-plane distance between two corresponding points is defined as the orthogonal distance from a point to the fitting plane of another point. For the point clouds to be stitched, estimate the transformation parameters by least squares adjustment, minimize the sum of the squares of the point-to-plane distances, and finally perform transformation to transform the point cloud using the estimated parameters for point cloud stitching.

2. The multi-view underwater three-dimensional point cloud stitching method based on sparse and dense point cloud fusion according to claim 1, characterized in that: The surface images of the underwater structure are obtained by collecting with the left camera and the right camera. Two cameras with different directions simultaneously capture the surface of an object in the same scene in different directions, and the morphology of the object surface is reconstructed through the pictures of the two cameras.

3. The multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds according to claim 2, characterized in that, Generate a dense point cloud based on stereo matching and triangulation, and output the correspondence between the dense point cloud and the pixels of the left camera image, specifically as follows: (31) Calibrate the binocular cameras using the checkerboard calibration method to obtain the internal and external parameters of the binocular cameras: (31.1) Camera calibration is the process of establishing the mutual position correspondence between the pixel coordinate system and the object points in space, and is also the process of determining the camera internal parameter matrix K, the external parameter matrix Rt, and the distortion coefficients. The internal parameter matrix K of the camera contains the basic attributes of the camera and is used to describe the optical characteristics of the camera. Its form is: Wherein: f x , f y , f s represents a scaling ratio relationship, which is calculated by the image distance f , the pixel size s and the true angle θ of the pixel arrangement in the camera plane: Wherein: s x , s y are the width and height of the pixel size respectively; (c x , c y ) represents the pixel coordinates of the camera principal point, c x is the pixel x-axis coordinate, c y is the pixel y-axis coordinate; f s From the image distance during imaging f , the pixel size in the horizontal direction s x and the true included angle θ of the pixel arrangement on the sensor plane are calculated: △r is the mathematical model of radial distortion, expressed as: Decompose the radial distortion into the x and y axis directions to obtain the Gaussian radial distortion model, where r is the distance from the image point to the center of the image plane, r 2 =x 2 +y 2 , k1, k2, and k3 represent the radial distortion coefficients, x and y respectively represent the abscissa and ordinate of the points in the ideal image coordinate system, and △X r represents the Gaussian radial distortion decomposed into the x-axis direction, and △Y r represents the Gaussian radial distortion decomposed into the y-axis direction, and are respectively expressed as: P(r) is the mathematical model of tangential distortion, where r is the distance from the image point to the center of the image plane, and r 2 = x 2 + y 2 , p1 and p2 represent the tangential distortion coefficients, expressed as: Decompose the tangential distortion into the x and y axis directions. x and y represent the abscissa and ordinate of points in the ideal image coordinate system respectively, and △x d represents the tangential distortion decomposed into the x-axis direction, and △y d represents the tangential distortion decomposed into the y-axis direction, and are respectively expressed as: In addition, the extrinsic parameters of the camera include a 3×3 rotation matrix R and a 3×1 translation matrix M to complete the affine transformation from the world coordinate system ( O w - x w y w z w ) to the camera coordinate system ( O - xyz ), where x w , y w , z w represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the point ( x w ,[[ID=2?]] y w , z w ) in the three-dimensional world coordinate system respectively, x , y , z represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the point ( x , y , z ) in the three-dimensional camera coordinate system respectively, and the affine transformation is expressed as follows: Matrix A is composed of a rotation matrix R and a translation matrix M, where R 11 , R 12 , R 13 , R 21 , R 22 , R 23 , R 31 , R 32 , R 33 respectively represent the corresponding position elements in the rotation matrix R, and M x , M y , M z respectively represent the corresponding position elements in the translation matrix M; And in the second step, the camera coordinate system needs to be converted to the image coordinate system ( O i – x i y i ), where x i , y i respectively represent the x-axis coordinate and y-axis coordinate of the point ( x i , y i ) in the two-dimensional image coordinate system, and the coordinate transformation is expressed as: In the formula, f is the image distance during imaging, α is the scale factor of the coordinates, and the B matrix is composed of the image distance f . Finally, the image coordinates are transformed into the pixel coordinate system ( O s – x s y s ), where x s , y s respectively represent the x-axis coordinate and y-axis coordinate of the point ( x s , y s ) in the pixel coordinate system: In the formula, s x , s y are the width and height of the pixel size; (c x , c y ) represents the pixel coordinates of the camera principal point, θ is the true angle of the pixel arrangement in the camera plane, the C matrix consists of 1 / s x , -cotθ / s x , c x , 1 / ( s y sinθ), c y . By converting and combining all coordinate systems and substituting them, the camera imaging mathematical model is obtained: ； (31.2) The internal parameter matrices K and distortion coefficients of the left and right cameras correct the original images, and the two corrected images lie in a plane and are parallel to each other; For a point (I, J, K) in the camera coordinate system, we have: where k1, k2, and k3 represent radial distortion coefficients, and p1, p2 represent tangential distortion coefficients, where R 11 ,R 12 ,R 13 ,R 21 ,R 22 ,R 23 ,R 31 ,R 32 ,R 33 represent the corresponding elements in the rotation matrix R respectively, M x ,M y ,M z represent the corresponding elements in the translation matrix M respectively, r is the distance from the image point to the center of the image plane, r 2 =x 2 +y 2 , f x , f y , f s ,c x and c y are elements in the internal parameter matrix K, x ’ and y ’ are the abscissa and ordinate of the camera coordinate point after z - coordinate normalization, x ” and y ” are the abscissa and ordinate of the distorted camera coordinate point, (x s ,y s ) are the pixel coordinates of the image point. The internal and external parameters are calibrated using a checkerboard calibration board. For each feature point on the calibration board, the following equation is satisfied: where r1, r2, and r3 are the three column vectors of the rotation matrix R, t represents the translation matrix M, K is the intrinsic parameter matrix, α is the scale factor of the coordinates, (x w , y w , z w ) is the coordinate in the world coordinate system of the point on the checkerboard, and establish the homography mapping relationship directly between the calibration board point (x w , y w , z w ) and the image point (x s , y s ): Among them, H is a homography matrix, and h 11 , h 12 , h 13 , h 21 , h 22 , h 23 , h 31 , h 32 respectively represent the corresponding position elements in the homography matrix H, s is an unknown scale factor, and H ’ = [h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 T is the matrix operation equation:​ (31.3) Use the least squares method to solve the homography matrix H. Since , where h1, h2, and h3 are the three column vectors of the homography matrix H, λ is the scale factor, r1 and r2 are the column vectors of the first two columns of the rotation matrix R, and t is the translation matrix T. Since the rotation matrix is an orthogonal matrix, its rows or column vectors are pairwise orthogonal and unit vectors. Thus, two constraint equations are obtained: K is the camera internal parameter matrix, h1 T , h2 T respectively represent the transposes of the column vectors h1 and h2, K -T represents the transpose of the inverse matrix of the K vector, K -1 represents the inverse matrix of the K vector. Let: B is a symmetric matrix, B 11 , B 12 , B 13 , B 21 , B 22 , B 23 , B 31 , B 32 , B 33 respectively represent the corresponding position elements in matrix B. Define the b matrix with a 6-dimensional vector. The T in the upper right corner of the matrix represents the transpose of the matrix: It can be seen that: h1 and h2 are the first two column vectors of the homography matrix H, and h i h1 and h2 can be taken according to the value of i, and h j h1 and h2 can be taken according to the value of j, and h i T The h represented i Inverse matrix, where v ij T The vectors are represented as follows. The T in the upper right corner of the matrix represents the transpose of the matrix: h 11 ,h 12 ,h 13 ,h 21 ,h 22 ,h 23 represent the corresponding position elements in the homography matrix H respectively, then the constraint equation is written in matrix operation form using matrix v ij and b: Among them, v 11 T represents the v ij T corresponding position element, v 12 , v 22 , respectively represent v ij corresponding position element, obtained from the matrix operation expression, multiple groups of v ij are substituted into the expression to be written in the form of the solution of the homogeneous equation: Vb = 0; the least squares solution of the homogeneous equation system is V T the eigenvector corresponding to the smallest eigenvalue of V, the eigenvector is normalized to obtain the required vector b, and then the symmetric matrix B is obtained. To obtain the symmetric matrix B, under the condition of a scale factor λ difference, the internal parameter is solved: Among them, B 11 , B 12 , B 13 , B 21 , B 22 , B 23 , B 31 , B 32 , B 33 respectively represent the elements at the corresponding positions in matrix B. f x , f y , f s , c x and c y are the elements in the intrinsic parameter matrix K. f x , f y , f s represents a scaling ratio relationship. (c x , c y ) represents the pixel coordinates of the camera principal point. After solving the intrinsic parameter matrix K, the extrinsic parameter matrix of each image is solved using the homography matrix H and the intrinsic parameter matrix K of each image: In the formula, λ is the scale factor, r1, r2, and r3 are the three column vectors of the rotation matrix R, h1, h2, and h3 are the three column vectors of the homography matrix H, t represents the translation matrix M, and K is the internal parameter matrix; Solving for distortion parameters: Let (x s1 , y s1 ) be the ideal image point coordinates without considering distortion, (x s , y s ) be the actual image point coordinates, x ’ and y ’ are the abscissa and ordinate of the camera coordinate point after z - coordinate normalization, x ” and y ” are the abscissa and ordinate of the camera coordinate point after distortion, f x , f y , f s , c x and c y are the elements in the internal parameter matrix K: Subtracting gives: Rearranging into the matrix operation form regarding distortion: Among them, the X matrix is the distortion parameter, k1, k2, and k3 represent the radial distortion coefficients, and p1, p2 represent the tangential distortion coefficients; the D matrix is known; in the E matrix, (x s1 , y s1 ) are the ideal image point coordinates, (x s , y s ) are the actual image point coordinates, and the T in the upper right corner of the matrix represents the transpose of the matrix. Therefore, the E matrix is also known; When there are multiple points, use the matrix least squares algorithm to solve for the camera distortion parameters k1, k2, k3, p1, and p2; (31.4) Use all the coefficients to correct the original images, and the two corrected images lie in a plane and are parallel to each other; (32) Generate a dense point cloud based on stereo matching and triangulation: By performing feature matching and disparity calculation on the images captured by the left and right cameras, obtain the disparity value of each pixel point in the scene; Utilize the disparity information, combined with the parameters of the camera, and calculate the coordinates of the object in three-dimensional space through triangulation; The two cameras on the left and right observe a point P in space simultaneously w with coordinates (x w , y w , z w ). O l and O r are the optical centers of the left and right cameras respectively. P l and P r are the image point pairs of the spatial point P w after epipolar alignment in the left and right cameras respectively. P l has coordinates (x l , y l ), and P r has coordinates (x r , y r ), where epipolar alignment is a method used in binocular stereo vision to increase the efficiency of searching for matching points; Let the distance between the optical centers of two cameras be \(b\), the focal lengths of the cameras be \(f\), and \(z\) be the depth information of a spatial point to the imaging plane. Then \(P\) w and \(P\) l 、 \(P\) r The relationship is obtained according to the similarity relationship of triangles: Let d = x l -x r We get: When obtaining the disparity d of a point in the two cameras of the binocular stereo vision system, the depth information of this point can be obtained, that is, the three-dimensional coordinates are restored; (33) Output the correspondence between the dense point cloud and the pixels of the left camera image; Using the left camera image as a template, the coordinates of the three-dimensional space points all have corresponding pixel coordinates in the image.

4. The multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds according to claim 1, wherein, Based on the received left camera image, use the structure from motion algorithm to generate a sparse point cloud using multi-view triangulation and output the correspondence between the sparse point cloud and the pixels of the left camera image, specifically as follows: (41) Detect and match feature points and perform matching and geometric verification. Use the scale-invariant feature transform algorithm to extract features, represent the feature points with 128-dimensional numerical descriptors, match the feature points, screen the matching point pairs based on the similarity of the descriptor vectors, and eliminate the abnormal matching point pairs through geometric verification; (42) Use multi-view triangulation to generate a sparse point cloud. Initialize and select two matching images, set the pose of one image as the identity matrix, estimate the transformation matrix using the matching point pairs, decompose to obtain the pose of the other image, and calculate the coordinates of the matching feature points in three-dimensional space through triangulation based on the camera poses of the two images to generate an initial sparse point cloud. As more images are added, match the new image with the feature points in the existing point cloud, and iteratively optimize the poses of all images and the three-dimensional point coordinates to gradually improve the sparse reconstruction of the scene; X W represents a three-dimensional point in space, O1 and O2 are the optical centers of two cameras, and the plane formed by these two camera optical centers and point X W is called the epipolar plane, and I1 and I2 are the physical image planes corresponding to these two cameras; In these two cameras, point X W has pixel coordinates x1 and x 2， respectively, which represent two projection observations of the 3D point X W ; x1' and x2' are the normalized plane coordinates corresponding to x1 and x2 respectively; the internal parameter matrices of the two cameras are K1 and K2; the coordinate transformation matrix of the camera coordinate system is R t, and there is a constraint: Among them, E is the essential matrix and F is the fundamental matrix. After completing feature extraction and matching, the symbols in the upper right corner of the matrix all represent matrix operations. Calculate the essential matrix and the fundamental matrix according to the above constraint conditions, and then obtain the external parameter information of the camera; Let point X W The pixel coordinates of W are x1 = (u1, v1, 1) and x2 = (u2, v2, 1), which are the feature points matched in two images. Here, u1, v1, u2, and v2 are the coordinates corresponding to the x-axis and y-axis respectively, and we get: Among them, F 11 , F 12 , F 13 , F 21 , F 22 , F 23 , F 31 , F 32 , F 33 respectively represent the elements at the corresponding positions in the fundamental matrix F. Let f = (F 11 , F 12 , F 13 , F 21 , F 22 , F 23 , F 31 , F 32 , F 33 ) to obtain: Each pair of matching points can provide a constraint for calculating the F matrix. Therefore, at least 8 pairs of matching feature points are required to solve for F. When the number of matching feature points is greater than 8, the least squares method or the RANSAC algorithm is used to solve for the F matrix. Given the camera parameters and the matching points, the coordinates of the 3D points are obtained through triangulation. Let the projection matrix of the i-th camera be P i =K i [R i ,t i , P i is the projection matrix of the i-th camera, K i is the internal parameter matrix of the i-th camera, R i is the rotation matrix of the i-th camera, t i is the translation matrix of the i-th camera. The coordinates of the 3D point in space are X W =[x,y,z,1] T , and the image coordinates of the corresponding feature point are x i =[u i ,v i ,1] T , u i , v i are the abscissa and ordinate of the x i point, and d i represents the depth of X W in the i-th view, that is, the coordinate z in the camera coordinate system c Then: Cross multiply both sides by x i We get: Further, we get: P i1 ,P i2 ,P i3 are respectively the horizontal vectors at the corresponding positions. The rank of the system of equations is 2, and at least one pair of matching feature points is required to solve the coordinates of the 3D point X; i ​ (43) The correspondence between the sparse point cloud and the pixels of the camera image is based on the input left camera image as a template, and the coordinates of the three-dimensional space points all have corresponding pixel coordinates in multiple images.

5. The multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds according to claim 4, characterized in that According to the obtained correspondences between the dense point cloud and the sparse point cloud and the pixels of the image respectively, calculate the correspondence between the dense point cloud and the sparse point cloud to obtain the initial alignment of the point cloud positions and complete the rough stitching of the point cloud, specifically as follows: Each three-dimensional point in the dense point cloud corresponds one-to-one with the pixel coordinates in the left camera image, and each three-dimensional point in the sparse point cloud corresponds one-to-one with the pixel coordinates in the input image. By using the two-dimensional coordinates of the image pixels as the index values, the three-dimensional points in the dense point cloud and the three-dimensional points in the sparse point cloud can be associated, so as to obtain the spatial position correspondence between the dense point cloud and the sparse point cloud; For each three-dimensional point X in the sparse point cloud xs , whose coordinates are (x xs , y xs , z xs ), there exists a pixel point X t in the left camera image corresponding to it. The coordinates of X t are (x t , y t ). Similarly, this pixel point X t must also correspond to a certain three-dimensional point X cm in the dense point cloud. The coordinates of X cm are (x cm , y cm , z cm ). Then, taking X t as the index, match the corresponding three-dimensional points in the sparse point cloud and the dense point cloud, calculate the difference between the point cloud coordinates, and take the average of the coordinate differences of the matching points to obtain the translation matrix M q . The values of x q , y q , and z q are the translation amounts corresponding to the x, y, and z directions respectively. Among them, n represents the total number of pixel points X t , and t is the index number of the pixel point: where x xst ,x cmt ,y xst ,y cmt ,z xst ,z cmt respectively represent the t-th element of x xs ,x cm ,y xs ,y cm ,z xs ,z cm ; Based on the spatial position correspondence between the dense point cloud and the sparse point cloud, a global translation matrix is calculated to align the dense point cloud as a whole to the position of the sparse point cloud, and the rough stitching of the point cloud is completed.

6. A multi-view underwater three-dimensional point cloud stitching system based on the fusion of sparse and dense point clouds, which is used to implement the multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds according to any one of claims 1 to 5, characterized in that, Including: A dense point cloud reconstruction module, configured to receive the surface images of an underwater structure binocular camera, where the binocular camera includes a left camera and a right camera, generate a dense point cloud based on stereo matching and triangulation, and output the correspondence between the dense point cloud and the pixels of the left camera image; A sparse point cloud reconstruction module, configured to generate a sparse point cloud based on the received left camera image by using the structure from motion algorithm and multi-view triangulation, and output the correspondence between the sparse point cloud and the pixels of the left camera image; A rough stitching module, configured to calculate the correspondence between the dense point cloud and the sparse point cloud according to the correspondences between the obtained dense point cloud and sparse point cloud and the image pixels respectively, obtain the preliminary alignment of the point cloud positions, and complete the rough stitching of the point cloud; A fine stitching module, based on the rough stitching of the point cloud, uses the iterative closest point algorithm to perform precise stitching on the point cloud and minimize the error in the overlapping area of the point cloud.

7. The multi-view underwater three-dimensional point cloud stitching system based on the fusion of sparse and dense point clouds according to claim 6, characterized in that: The dense point cloud reconstruction module includes a binocular stereo vision system, specifically set as a binocular camera, which includes two cameras with different directions, simultaneously takes pictures of the surface of an object in the same scene in different directions, and reconstructs the morphology of the object surface through the pictures of the two cameras.

8. A terminal device, comprising a memory, a processor, and a computer program stored in the memory and capable of running on the processor, characterized in that, The memory stores a computer program that can run on a processor. When the processor loads and executes the computer program, the multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds described in any one of claims 1 to 5 is adopted.

9. A storage medium containing computer-executable instructions, characterized in that, The computer-executable instructions are used to execute the multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds described in any one of claims 1 to 5 when executed by a computer processor.

Images