Multi-view underwater three-dimensional point cloud splicing method based on sparse and dense point cloud fusion

Abstract

The invention discloses a multi-view underwater three-dimensional point cloud splicing method based on sparse and dense point cloud fusion, and relates to the technical field of structure underwater detection. The method comprises the steps that a binocular camera surface image of an underwater structure is received, the binocular camera comprises a left camera and a right camera, dense point clouds are generated based on stereo matching and triangulation, and the corresponding relation between the dense point clouds and image pixels of the left camera is output; and based on the received left camera image, a sparse point cloud is generated by adopting a motion structure recovery algorithm and utilizing multi-view triangulation, and a corresponding relation between the sparse point cloud and the pixels of the left camera image is output. According to the method, when the point cloud of multi-view reconstruction is spliced, the surface features of the underwater structure can be well restored, feature extraction is not needed, the method has high accuracy and robustness, and high-precision point cloud splicing can be achieved in the underwater complex environment, so that the method has good application prospects in actual underwater structure detection and image splicing 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, which 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 limited by underwater lighting and imaging quality, it is difficult to obtain stable global features. 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 object, 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, including the following steps: Receiving the 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; 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 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, using the Iterative Closest Point algorithm to perform precise stitching on the point cloud to minimize the error in the overlapping area of the point cloud.

[0005] Further, the surface images of the underwater structure are collected by the left camera and the right camera. Two cameras with different directions simultaneously take pictures of the surface of an object in the same scene from different directions, and the morphology of the object surface is reconstructed through the pictures of the two cameras.

[0006] Furthermore, a dense point cloud is generated based on stereo matching and triangulation, and the correspondence between the dense point cloud and the pixels of the left camera image is output as follows: (31) Calibrate the binocular camera using the checkerboard calibration method to obtain the internal and external parameters of the binocular camera: (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 it 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: Where: f x , f y represents a scaling ratio relationship, which is calculated through the image distance f , the pixel size s , and the true included angle θ of the pixel arrangement on the camera plane: Where: 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, and c y is the pixel y-axis coordinate; f s is calculated from the image distance f , the horizontal pixel size s x , and the true included angle θ of the pixel arrangement on the sensor plane: △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. 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 points in the ideal image coordinate system, △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: 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, which are expressed as: The tangential distortion is decomposed into the x and y axis directions. x and y respectively represent the abscissa and ordinate of the points 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 they are respectively expressed as: 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: 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, 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, the B matrix is composed of the image distance f , and 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: 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, the camera imaging mathematical model is obtained: ; (31.2) The internal parameter matrix K and distortion coefficients of the left and right cameras correct the original images, and the two corrected images are in a plane and parallel to each other; For a point (I, J, K) in the camera coordinate system, there is: Among them, k1, k2, and k3 represent radial distortion coefficients, and p1 and p2 represent tangential distortion coefficients. 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 , 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 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 ) is the pixel coordinate of the image point. The internal and external parameter values are calibrated by using a checkerboard calibration board. For each feature point on the calibration board, the following formula is satisfied: In the formula, 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, α is the scale factor of the coordinates, (x w , y w , z w) is the coordinate in the world coordinate system of the points on the chessboard, and establish the homography mapping relationship between the calibration board points (x w , y w , z w ) and the image points (x s , y s ): where, 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 matrix operation equation of : (31.3) Use the least squares method to solve the homography matrix H, and because , where h1, h2, h3 are the 3 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, 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, and 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 elements in matrix B, define matrix b with a 6-dimensional vector, and 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, and 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 respectively represent the corresponding elements in the homography matrix H, then the constraint equation is written in matrix operation form with matrix v ij and b: Among them, v 11 T represents the corresponding element of v ij T The corresponding element, v 12 , v 22 , respectively represent the corresponding 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 The eigenvector corresponding to the smallest eigenvalue of V, normalize the eigenvector to obtain the required vector b, and then obtain matrix B. Obtain matrix B. Under the condition of a scale factor λ, solve the internal parameter: Among them, 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, fx , f y , f s , c x and c y are 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 is solved using the homography matrix H and the internal parameter matrix K of each image: In the formula, λ 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; 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, and 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 elements in the internal parameter matrix K: Subtract to get: Arrange it 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 ) is the ideal image point coordinate, (x s , y s ) is the actual image point coordinate, 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, the values of the camera distortion parameters k1, k2, k3, p1, and p2 are solved using the matrix least squares algorithm; (31.4) Correct the original images using all the coefficients. The two corrected images are located 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 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 the parameters of the camera, the coordinates of the object in three-dimensional space are calculated through triangulation; (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 respectively. The coordinate of P l is (x l , y l ), and the coordinate of P r is (x r , y r ). Among them, epipolar alignment is a method used in binocular stereo vision to increase the search efficiency of matching points; 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 triangle similarity relationship: Let d = x l - x r to get: 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; (33)Output the correspondence between the dense point cloud and the pixels of the left camera image; Taking the left camera image as a template, the coordinates of the three-dimensional points in the three-dimensional space all have corresponding pixel coordinates in the image.

[0007] Furthermore, based on the received left camera image, the Structure from Motion (SfM) 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: (41)Detect and match feature points, perform matching and geometric verification, extract features using the Scale-Invariant Feature Transform (SIFT) algorithm, 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)Generate a sparse point cloud by multi-view triangulation. 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 by 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 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 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; 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 as follows: Among them, E is the essential matrix, F is the fundamental matrix. After feature extraction and matching, the symbols in the upper right corner of the matrix all represent matrix operations. According to the above constraint conditions, calculate the essential matrix and the fundamental matrix, and then obtain the external parameter information of the camera; Let point X WThe pixel coordinates x1 = (u1, v1, 1) and x2 = (u2, v2, 1) are the feature points matched in two images, where u1, v1, u2, and v2 are the coordinates corresponding to the x-axis and y-axis respectively, and we have: Among them, F 11 ,F 12 ,F 13 ,F 21 ,F 22 ,F 23 ,F 31 ,F 32 ,F 33 respectively represent the corresponding 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: 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 , where 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 c in the camera coordinate system. Then: Cross multiply both sides by x i We get: Furthermore, we get: 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; (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.

[0008] 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, and obtain the preliminary alignment of the point cloud positions to complete the rough stitching of the point cloud, as follows: (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 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 xs in the sparse point cloud, its coordinates are (x xs , y xs , z xs ). 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 ). 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, using 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 find the average value of 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. 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 and x cmt , y xst , and y cmt , z xst , and z cmt represent the t-th elements of x xs , x cm , y xs , and y cm , z xs , and z cm respectively; (52) 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, completing the rough stitching of the point cloud.

[0009] Further, on the basis of the rough stitching of the point cloud, the iterative closest point algorithm is used to accurately stitch the point cloud, minimizing the error in the overlapping area of the point cloud, specifically as follows: (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. Then perform matching 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; 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: For point X v with coordinates (x v , y v , z v ), its voxel index calculation is as follows: Traverse all voxels. If there are points from both P1 and P2 in the voxel, then the voxel belongs to the overlapping area V overlap : In V overlap perform uniform sampling to obtain a sparse but evenly distributed point set S. Let P be the point set of the fixed point cloud. For each voxel vi，j，k ∈V overlap , uniformly select a representative point X within this voxel s Form point set S: 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 another point cloud. For each sampling point X ’ in P s , find the nearest neighbor point X ’ 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 the threshold τ d for rejection: If d s > τ d , then consider the match incorrect and reject this point pair; 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: Regard 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 P to be registered = {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 vector N jThe corresponding coordinate values, where 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 in the target point cloud Q j , and obtain the normal vector N j . For each matching point pair (X i , Y j ), calculate the point-plane distance error d i : where N j T represents the transpose matrix of N j . Parameter optimization is performed by minimizing the weighted sum of squares of point-plane distances for all corresponding points: 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: where is the median and σ k is the noise estimate for the k-th point cloud pair; Let the parameter vector be β. 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 residual distance, linearized into a f function related to x0: where d i (x) is the point-plane distance error equation with respect to x, vi represents the adjusted residual 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 residual distance, and △β is the parameter correction; Solve for the parameter correction △β iteratively: In the formula, β0 is the initial parameter value, A is the design matrix, and the weight matrix P = diag(ω1, ω2…ω n ), ω i is the corresponding weight; The iterative update formula is: Converge when |△β| is below 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 mis - matched points. The internal iteration includes linearizing the least - squares solution for △β, 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, the transformation parameters are estimated through least - squares adjustment to minimize the sum of the squares of the point - to - plane distances. Finally, a transformation is performed, and the point clouds are transformed using the estimated parameters for point cloud stitching.

[0010] The present invention has at least the following beneficial effects: 1. When the present invention stitches the point clouds for multi - view reconstruction, it can better restore the surface features of the underwater structure without feature extraction, with 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.

[0011] 2. The present invention does not need to rely on feature markers and can complete the full - view reconstruction of complex underwater structures in a single acquisition only through binocular vision devices.

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

[0013] Of course, any product implementing the present invention does not necessarily need to achieve all the above - mentioned advantages simultaneously. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 is a flowchart of the method of the present invention; Figure 2 is a flowchart of dense point cloud reconstruction in the present invention; Figure 3 is a schematic diagram of the binocular camera vision principle in the present invention; Figure 4 is a flowchart of the three - dimensional reconstruction of the structure - from - motion algorithm in the present invention; Figure 5 is a flowchart of rough stitching based on pixel index values in the present invention; Figure 6 is a diagram of the refined stitching result of multi - view underwater three - dimensional point clouds in the present invention; Figure 7 is a schematic diagram of the fast stitching result of the underwater structure surface 3D reconstruction point cloud in the present invention; Figure 8 is a schematic diagram of the structure of the system of the present invention; Figure 9 is a schematic diagram of the implementation process in the embodiment of the present invention. Detailed implementation manners

[0015] The following will clearly and completely describe the technical solutions in the embodiments of the present disclosure with reference to the accompanying drawings in the embodiments of the present disclosure. Obviously, the described embodiments are only a part rather than all of the embodiments of the present disclosure. 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.

[0016] Embodiment 1: 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.

[0017] 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: S1. Receive the surface images of the 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: (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; (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; 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: Where: f x , f y represents a scaling ratio relationship, which is calculated through the image distance f , pixel size s and the true included angle θ of the pixel arrangement on the camera plane: Where: s x , s y is the width and height of the pixel size; (c x , c y ) represents the pixel coordinates of the camera principal point; 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, calculate: △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. 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: 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: Decompose the tangential distortion 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 are respectively expressed as: In addition, the external 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 respectively represent the point in the three-dimensional world coordinate system ( xw , y w , z w ), the x-axis coordinate, y-axis coordinate, z-axis coordinate, x , y , z respectively represent the x-axis coordinate, y-axis coordinate, z-axis coordinate of the point ( x , y , z ) in the three-dimensional camera coordinate system. The affine transformation is expressed as follows: 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; 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, y-axis coordinate of the point ( x i , y i ) in the two-dimensional image coordinate system. 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 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: 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, and the C matrix is composed of 1 / s x , -cotθ / s y , 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: (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 in a plane and parallel to each other; For the camera coordinate system point (I, J, K), there is: Among them, k1, k2, and k3 represent the radial distortion coefficients, p1 and p2 represent the 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 , fx , 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 by relying on the checkerboard calibration board. For each feature point on the calibration board, the following formula is satisfied: In the above formula, r1, r2, r3 are the three column vectors of the rotation matrix R, t represents the translation matrix M, K is the internal parameter matrix, α is the scale factor of the coordinate, (x w , y w , z w ) are the coordinates 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 ): 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 coefficient, and H ’ = [h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 T The matrix operation equation of: As can be seen from the above, the equations of multiple corresponding points are superimposed and can be 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 are required in a single photo to solve for the homography matrix H. There are 8 equations to solve for 8 unknowns, and each image has a homography matrix; And 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 are unit vectors. Thus, two constraint equations are obtained: 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: 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 6 - dimensional vector. The T in the upper right corner of the matrix represents the transpose of the matrix: It can be known that: 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: h 11 h 12 h 13 h 21 h 22 h 23Denote the corresponding position elements in the homography matrix H respectively, then the constraint equation is written in matrix operation form with matrix v ij and b: where, v 11 T denotes the corresponding position element of v ij T ; v 12 , v 22 respectively denote the corresponding position elements of v ij . It can be seen from the above formula that 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. Normalize the eigenvector to obtain the required vector b, and then obtain matrix B. There are 6 unknowns for vector b, and at least 3 groups of v ij are required, that is, 3 H matrices, and 6 equations are used to solve 6 unknowns; therefore, at least 3 different calibrated pictures are required to solve matrix b. After obtaining matrix b, construct matrix B. Under the condition of a scale factor difference (B = λK -T K -1 , λ is the scale factor), solve the internal parameter: where B 11 , B 12 , B 13 , B 21 , B 22 , B 23 , B 31 , B 32 , B 33 respectively denote the corresponding position 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: where λ 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 parameter matrix; Next, solve for 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 point after distortion, f x , f y , f s , c x and c y are the elements in the intrinsic parameter matrix K, then: Subtracting gives: Rearranging into the matrix operation form regarding distortion: where the X matrix is the distortion parameter to be solved, k1, k2, and k3 represent the radial distortion coefficients, and p1, p2 represent the tangential distortion coefficients; The camera intrinsic and extrinsic parameters required in the D matrix have been obtained previously, 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 arranged into the matrix operation form, and the values of the camera distortion parameters k1, k2, k3, p1, and p2 can be solved using the matrix least - squares algorithm; Finally, use all the coefficients to correct the original image, and the two corrected images lie in a plane and are parallel to each other; (S12) While ensuring the overlap of the image fields of view, 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: 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; Using the disparity information and combining it with the parameters of the camera, the coordinates of the object in three-dimensional space are calculated through triangulation; The left and right cameras simultaneously observe a point P in space 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. Epipolar alignment is a method used in binocular stereo vision to increase the search efficiency of matching points. Assuming the distance between the optical centers of the two cameras is b and the focal lengths of the cameras are both f, then the relationship between P w and P l , P r can be obtained according to the similarity relationship of triangles: Let d = x l - x r It can be obtained: When the disparity d of a certain point in the binocular stereo vision system's two cameras is obtained, the depth information of this point can be obtained, that is, the three-dimensional coordinates are restored; (S13) Output the correspondence between the dense point cloud and the pixels of the left camera image; Taking the left camera image as a template, the coordinates of the three-dimensional space points all have corresponding pixel coordinates in the image; Specifically, the specific acquisition device for the underwater stereo vision image is an underwater binocular camera, which uses the optical method of the camera 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; Regarding the technical solution of this embodiment, such as Figure 2The following is a flowchart for dense point cloud reconstruction. It mainly utilizes the binocular vision principle and two images with overlapping regions collected by a binocular camera to reconstruct the surface topography of an underwater structure. The main processes include underwater calibration of the calibration board, selection of the image overlapping region, and calculation of the coordinates of the dense point cloud. For the binocular system calibration using the calibration board, a checkerboard calibration board is used. Both the binocular system and the calibration board are placed in an underwater environment. Calibration images of the calibration board are collected at different object distances to calculate the internal and external parameters of the left and right cameras. The epipolar correction and distortion correction are performed using the internal and external parameters. 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 parameters of the camera, 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; For 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 ), the imaging point P1, and the actual point P are collinear. Similarly, the optical center of camera B (O c2 ), the imaging point P2, and the 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; S2. Based on the received left camera image, use the structure from motion recovery algorithm 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: (S21) 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; (S22) Generate a sparse point cloud using multi-view triangulation. Initialize by selecting two matching images, setting the pose of one image to the identity matrix, estimating the transformation matrix using the matching point pairs, decomposing to obtain the pose of the other image, and calculating the coordinates of the matching feature points in 3D 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 images 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; X W represents a 3D point in space, and O1 and O2 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; I1 and I2 are the physical image planes corresponding to these two cameras. 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 3D 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: 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 the fundamental matrix can be calculated according to the above constraint conditions, and then the external parameter information of the camera can be obtained; Let the pixel coordinates of point X W be x1 = (u1, v1, 1), x2 = (u2, v2, 1) as the matching feature points of two images, where u1, v1, u2, v2 are the coordinate values corresponding to the x-axis and y-axis respectively, and we get: 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 , F 13 , F 21 , F 22 , F 23 , F 31 , F 32 , F 33) It follows that: 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 there are more than 8 matching feature points, the least squares method or the RANSAC algorithm is used to solve the F matrix. Given the known camera parameters and 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 The corresponding image coordinates of the feature point are x i =[u i , v i , 1] T where 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 It can be obtained that: Furthermore, it can be obtained that: P i1 , P i2 , P i3 are the horizontal vectors at the corresponding positions of P i The rank of this system of equations is 2. At least one pair of matching feature points is required to solve the coordinates of the 3D point X. In fact, there are multiple feature matching point pairs for a 3D point in the input images of multiple views. The coordinates of the 3D point X can be accurately solved by the least squares method or the RANSAC method; (S23) 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 3D space points have corresponding pixel coordinates in multiple images; Specifically, as Figure 4The following is 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 retrieval and matching, use the SIFT feature extraction method to find sparse feature points in the images, and describe the feature points using a numerical descriptor (128-dimensional vector). Then, match the image pairs again. Next, match the feature points on the matched images, and use the similarity measure of the 128-dimensional descriptor vectors of the feature points to perform feature point matching. Geometric verification is performed on 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 matching 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, and 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; S3. According to the correspondence 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, as follows: (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. Using the two-dimensional coordinates of the image pixels as the index values, the 3D points in the dense point cloud and the 3D 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 3D point X in the sparse point cloud xs (x xs , y xs , z xs ), there exists a pixel point X t (x t , y t ) in the left camera image corresponding to it. Similarly, this pixel point X t must also correspond to a certain 3D point X cm (x cm , y cm , z cm ) in the dense point cloud. Then, using X t as the index, match the corresponding 3D points in the sparse point cloud and the dense point cloud, calculate the difference between the point cloud coordinates, and calculate the average value of the coordinate differences of the matched 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 points: 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 ; (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; Specifically, as Figure 5 shown in the flowchart of rough stitching based on pixel index values, in the initial state, there are phenomena of coordinate system mismatch and too far distance between different reconstructed dense point clouds and SFM sparse point clouds. 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 index values, 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; 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, minimizing the error in the overlapping area of the point cloud, specifically as follows: (S41)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 to make the points in the overlapping area evenly distributed. Perform matching 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 ,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: The point X v (x v ,y v ,z v) The 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 : In V overlap perform uniform sampling 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 in this voxel to form the point set S: That is, select the centroid of all points in 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 ’ , find the nearest neighbor point X s in P ’ : ’ s : To improve the accuracy of the match, 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 , then consider the match incorrect and reject this point pair; (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: 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 , aim to find R icp and ticp , such that the transformed point cloud P ’ Optimal match Q: Estimate the normal vector N of the target point cloud Q j : 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, and 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 in the target point cloud Q j , and obtain the normal vector N j , for each matching point pair (X i , Y j ), calculate the point-plane distance error d i : where N j T represents the transpose matrix of N j , and parameter optimization is performed by minimizing the weighted sum of squares of the point-plane distances of all corresponding points: 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: where is the median, and σ k is the noise estimate of the kth point cloud pair; 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 residual distance, linearized into a f function related to x0: where 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 residual distance, and △β is the parameter correction; Solve for the parameter correction △β by iteration: where β0 is the initial value of the parameter, A is the design matrix, and the weight matrix P = diag(ω1, ω2…ω n ), ω i is the corresponding weight; The iterative update formula is: Until convergence, that is, |△β| 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 of △β, 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 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 to transform the point clouds using the estimated parameters for point cloud stitching; 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; Specifically, as Figure 6 is the fine - stitching result diagram of multi - view underwater three - dimensional point clouds: Initial state, the point clouds are already roughly aligned. A voxel hull represented by the low - resolution of the volume occupied by the point clouds is used to find and determine the overlapping area (searching for 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 within the overlapping area are evenly distributed. Then, matching is performed again to find the nearest points of the selected subset in another point cloud (matching of selected points and nearest points), and incorrect correspondences are rejected based on the compatibility of the points; 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, which will minimize the sum of the squares of the point - to - plane distances, and finally, a transformation is performed to transform the point clouds using the estimated parameters (minimization of the point - to - plane distance); The point clouds are stitched by transforming the point clouds using the estimated parameters. Finally, the convergence criterion is tested. If not satisfied, the process will start over until the result converges, and the optimal rotation and translation matrices are exported (fine - stitching result).

[0018] Next, the present invention will be further elaborated in combination with specific implementation cases: 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.

[0019] Figure 9 This is a 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 specifications 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.

[0020] Embodiment Two: 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 for implementing the above - mentioned multi - view underwater three - dimensional point cloud stitching method based on the fusion of sparse and dense point clouds, including: 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; 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; A rough stitching module, which is used to calculate the corresponding relationship between the dense point cloud and the sparse point cloud according to the corresponding relationships 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; The precise stitching module, based on the rough stitching of the point cloud, uses the iterative closest point algorithm to precisely stitch the point cloud and minimize the error in the overlapping area of the point cloud.

[0021] Embodiment 3: 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.

[0022] Embodiment 4: The present invention provides a storage medium containing computer-executable instructions. When the computer-executable instructions are executed by a computer processor, they 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.

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 on 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 using 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.

2. 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: The surface images of the underwater structure are acquired by the left camera and the right camera. Two cameras with different directions simultaneously capture the surface of an object in the same scene from 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 sparse and dense point cloud fusion according to claim 2, wherein, 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 it is also the process of determining the camera internal parameter matrix K, external parameter matrix Rt, and 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 included 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 and is 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, r 2 = x 2 + y 2 , and p1, p2 represent the tangential distortion coefficients, expressed as: Decompose the tangential distortion into the x and y axis directions. x and y respectively represent the abscissa and ordinate of a point in the ideal image coordinate system, 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 , 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 represent the corresponding position elements in the rotation matrix R respectively, and M x , M y , M z represent the corresponding position elements in the translation matrix M respectively; 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 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 to form. 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: 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 substituting all coordinate systems, the camera imaging mathematical model is obtained: ； (31.2) Use the internal parameter matrix K and distortion coefficients of the left camera and the right camera to correct the original images. The two corrected images are located in a plane and are parallel to each other; For the camera coordinate point (I, J, K), there is: 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 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: where 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 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 the 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 The matrix operation equation of:​ (31.3) The homography matrix H is solved 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 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 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: 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. 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 known that: h1 and h2 are the first two column vectors of the homography matrix H, h i h1 and h2 can be taken according to the value of i, h j h1 and h2 can be taken according to the value of j, 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 respectively represent the corresponding position elements in the homography matrix H, then the constraint equation is written in matrix operation form using the matrix v ij and b: Among them, v 11 T represents the v ij T corresponding position element, v 12 , v 22 , respectively represent the 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 minimum eigenvalue of V. Normalize the eigenvector to obtain the required vector b, and then obtain the symmetric matrix B. Obtain the symmetric matrix B. Solve the internal parameter under the condition of a scale factor λ difference: Among which 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, f x , f y , f s ,c x and c y are 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, and (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, and 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: Subtract to get: Rearrange it 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 the values of the camera distortion parameters k1, k2, k3, p1, and p2; (31.4) Use all the coefficients to correct the original images. The two corrected images are located in a plane and are parallel to each other; (32) Generate a dense point cloud based on stereo matching and triangulation: Through feature matching and disparity calculation of the images captured by the left and right cameras, obtain the disparity value of each pixel point in the scene; Use the disparity information and combine it with the parameters of the camera to calculate the coordinates of the object in the 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 ). Here, 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 the spatial point to the imaging plane. Then \(P\) w and \(P\) l 、 \(P\) r The relationship of \( Let d = x l -x r We get: When obtaining the disparity d of a certain 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; Taking 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, characterized in that Based on the received left camera image, a structure from motion recovery algorithm is used to generate a sparse point cloud using multi-view triangulation, and the correspondence between the sparse point cloud and the pixels of the left camera image is output 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 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; (42) Generate a sparse point cloud using multi-view triangulation. 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 the 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， 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 intrinsic 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 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; Set point X W The pixel coordinates of are x1 = (u1, v1, 1) and x2 = (u2, v2, 1), which are the feature points matched by two images. Among them, u1, v1, u2, and v2 are the coordinate values corresponding to the x-axis and y-axis respectively, and we get: Among them, F 11 and F 12 and F 13 and F 21 and F 22 and F 23 and F 31 and F 32 and 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, and 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 from the i-th perspective, that is, the coordinate z in the camera coordinate system c Then: Cross multiply both sides by x i We get: Furthermore: 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 for the coordinates of the 3D point X; i ​ (43) The correspondence between the sparse point cloud and the pixels of the camera image 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.

5. The multi-view underwater three-dimensional point cloud stitching method based on sparse and dense point cloud fusion according to claim 4, wherein, 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: (51) 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. Using the two-dimensional coordinates of the image pixels as 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, 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 , where x q , y q , z q are the translation amounts corresponding to the x, y, and z directions respectively. Here, 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 elements of x xs ,x cm ,y xs ,y cm ,z xs ,z cm ; (52) 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.

6. The multi-view underwater three-dimensional point cloud stitching method based on the fusion of sparse and dense point clouds according to claim 5, characterized in that 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 as follows: (61) 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 to find the closest points of the selected subset of points to be matched in another point cloud, and reject the 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 as follows: 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 at the same time, 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 , search for the nearest neighbor point X ’ in P ’ s : Use the geometric compatibility of points to filter out incorrect matches. Let d s be the distance between a pair of matching points (X s , X ’ S ). Set a threshold τ d for rejection: If d s > τ d , it is considered that there is a matching error and this point pair is removed; (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 precise stitching: Take the sparse point cloud as the target point cloud Q = {Y j = (x j , y j , z j ) T}, and the dense point cloud as the point cloud P to be registered = {X i = (x i , y i , z i ) T . The initial transformation parameters are the rotation matrix R o and the 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 . 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 with respect to x0 into a function related to f x0: 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 as follows: It converges 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 - plane distances, and filtering out mis - matched points. The internal iteration includes linearized least - squares solution of △β, updating 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 fitted plane of another point. For the point clouds to be spliced, the transformation parameters are estimated by least - squares adjustment, minimizing the sum of the squares of the point - to - plane distances. Finally, a transformation is performed, and the point clouds are transformed using the estimated parameters for point cloud splicing.

7. 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 6, characterized in that, It includes: A dense point cloud reconstruction module, which is used to receive the surface images of an underwater structure binocular camera. The binocular camera includes a left camera and a right camera. Based on stereo matching and triangulation, a dense point cloud is generated, and the correspondence between the dense point cloud and the pixels of the left - camera image is output; A sparse point cloud reconstruction module, which is used to generate a sparse point cloud based on the received left - camera image using a structure - from - motion algorithm with multi - view triangulation, and output the correspondence between the sparse point cloud and the pixels of the left - camera image; A rough splicing module, which is used 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 a preliminary alignment of the point cloud positions, and complete the rough splicing of the point cloud; A fine splicing module, which, based on the rough splicing of the point cloud, uses the iterative closest point algorithm to accurately splice the point cloud, minimizing the error in the overlapping area of the point cloud.

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

9. 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, it adopts the multi - view underwater three - dimensional point cloud splicing method based on the fusion of sparse and dense point clouds described in any one of claims 1 to 6.

10. 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 splicing method based on the fusion of sparse and dense point clouds described in any one of claims 1 to 6 when executed by a computer processor.

Images