Tunnel fracture identification method based on three-dimensional live-action reconstruction and orderly reacquisition of virtual camera

Abstract

The invention discloses a tunnel fracture recognition method based on three-dimensional live-action reconstruction and ordered re-acquisition of a virtual camera. The method specifically comprises the following steps: acquiring a multi-view rock fracture image; acquiring three-dimensional point cloud data of a tunnel wall surface according to the acquired rock mass fracture image, reconstructing a geometrical shape of the tunnel, mapping acquired texture information to a three-dimensional grid vertex through an image processing algorithm, and generating a three-dimensional tunnel model; setting parameters of the virtual camera, orderly collecting to form a virtual image, and enabling the position of the virtual camera to correspond to the real position of the model; and extracting information of the crack of the virtual image, and projecting the coordinates of the two-dimensional image back to the three-dimensional space through the projection matrix to obtain the space coordinates of the crack. A virtual image is generated by combining a texture mapping technology of a three-dimensional model, so that the limitation of a visual angle, illumination and space in actual image acquisition is overcome, and high-precision identification of tunnel cracks is realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of surrounding rock structure detection, and in particular relates to a tunnel fissure identification method based on three-dimensional real scene reconstruction and orderly re-collection by a virtual camera. Background Art

[0002] During tunnel construction, the underground environment is complex and ever-changing. The identification and recording of rock structure elements (such as cracks, joints, and bedding) in tunnel projects is crucial for tunnel design, construction, and maintenance. Traditional rock structure identification methods rely primarily on manual surveys, which suffer from low identification efficiency, irregular data recording, and subjective operator influence. Non-contact measurement technologies (3D laser scanning and photogrammetry) are increasingly being applied to fracture acquisition, but they still have significant limitations: high laser scanning equipment costs, low data processing efficiency, and sparse texture information. While photogrammetry can acquire high-resolution images, interpreting fracture spatial location is cumbersome, requiring multiple pre-calibration of the measurement instrument and the camera position, which affects fracture spatial location efficiency. Furthermore, due to limited shooting conditions, single-view photogrammetry often fails to provide a complete perspective and high-quality texture information. This limitation results in incomplete or inaccurate capture of tunnel surface features, which in turn affects fracture identification accuracy. It is proposed to first construct a three-dimensional real-scene model of the tunnel through on-site multi-view photography images, and then use virtual camera technology to simulate different perspectives and shooting conditions for orderly re-collection to obtain ideal images to improve recognition accuracy. At the same time, the known spatial position of the virtual camera facilitates the spatial coordinate positioning of the crack, thereby improving the efficiency of tunnel crack space positioning. Summary of the Invention

[0003] The purpose of the present invention is to provide a tunnel fissure identification method based on three-dimensional real scene reconstruction and orderly re-collection of virtual cameras, which can effectively improve the recognition accuracy and efficiency of rock mass fissures.

[0004] The technical solution adopted by the present invention is a tunnel fissure identification method based on 3D real scene reconstruction and orderly re-collection by a virtual camera, which is specifically implemented according to the following steps:

[0005] Step 1: Use the terminal to obtain multi-view rock fracture images;

[0006] Step 2: Based on the collected rock fracture images, obtain 3D point cloud data of the tunnel wall, reconstruct the tunnel geometry, and map the acquired texture information to 3D mesh vertices using an image processing algorithm to generate a 3D tunnel model.

[0007] Step 3: Set the parameters of the virtual camera, collect data in an orderly manner to form a virtual image, and match the position of the virtual camera with the real position of the model;

[0008] Step 4: Extract the information of the cracks in the virtual image and project the two-dimensional image coordinates back into the three-dimensional space through the projection matrix to obtain the spatial coordinates of the cracks.

[0009] The present invention is also characterized in that:

[0010] In step 2, specifically:

[0011] Step 2.1: For each rock fracture image, use the SIFT method to extract a set of key points with scale and rotation invariance, and generate a corresponding description vector for each key point, which is the feature point descriptor;

[0012] Step 2.2: For the descriptor sets in the rock mass fracture images A and B in adjacent frames, the feature point matching relationship between the images is calculated based on the similarity between the description vectors, and a 2D-2D corresponding point set between the images is established;

[0013] Step 2.3, use the RANSAC algorithm to estimate the essential matrix E1;

[0014] Step 2.4, perform singular value decomposition on the estimated essential matrix E1 to obtain the estimated essential matrix E2 = UDV T , where U and V are orthogonal matrices, D is a diagonal matrix, and the third singular value in D is set to zero, constructing the modified matrix D′=diag(s,s,0) and reconstructing the essential matrix E′=UD′V T At the same time, keep the U and V matrices;

[0015] Step 2.5, according to the formula Calculate the epipolar geometric error of all point pairs. When the error is less than the threshold of 1.5, the point is marked as an "inlier point", otherwise it is an "outlier point". After several rounds of iteration, the E2 with the most inliers is retained as the optimal solution. Output the estimated essential matrix E2 and the corresponding inlier point set.

[0016] Step 2.6: Based on the modified essential matrix E′=UD′V obtained in step 2.4 T , use the U and V matrices to decompose the relative posture of the camera, and construct four sets of candidate (R, t) combinations according to the standard decomposition method. Each set of rotation and translation corresponds to a different camera posture. The rotation selection matrix Rz represents the standard orthogonal matrix of ±π / 2 rotation around the z axis;

[0017] R1=UR Z (π / 2)V T ,t1=U[:,2]

[0018] R2=UR Z (-π / 2)V T ,t2=U[:,2]

[0019] R3=URZ (π / 2)V T ,t3=-U[:,2]

[0020] R4=UR Z (-π / 2)V T ,t4=-U[:,2]

[0021] In step 2.7, only one of the four solutions is correct, that is, under this combination, all the 3D points obtained by triangulation of the matching points are simultaneously in front of both cameras. The correct solution is determined by checking the 3D points, thus obtaining the pose matrix P = [R|t] for the current frame.

[0022] Step 2.8: Calculate the coordinates of the 3D points using triangulation. First, calculate the projection matrix M of each rock fracture image based on the known intrinsic and extrinsic matrix K of the camera. Let u be the homogeneous coordinates of the 2D point, U be the homogeneous coordinates of the 3D point, and M be the projection matrix. Then, the projection matrix satisfies u = MU.

[0023] Step 2.9, based on the projection matrix and corresponding points of the rock fracture images of adjacent frames, the projection matrix of the i-th camera is M i =K i [R i ,t i ]=[M i1 ,M i2 ,M i3 ] T , the corresponding pixel coordinate is x i =[u i ,v i ,1] T According to the projection equation, we can get d i x i =M i X, cross-multiply both sides of the equation by x i You can get x i ×M i X=0; rearranging the equations yields the scene's three-dimensional point cloud X;

[0024]

[0025] Step 2.10: Verify the four candidate poses according to the positive depth constraint, and finally select the one that satisfies all points with positive Z values ​​as the camera pose [R|t] of the current frame image.

[0026] Step 2.11: After obtaining preliminary 3D points through triangulation, perform bundle adjustment optimization to minimize the reprojection error of all 3D points and camera parameters.

[0027] Step 2.12: Based on the acquired sparse 3D point cloud X, the camera pose matrix [R|t], and the intrinsic parameter matrix K, we further calculate the projection matrix Mi = K[Ri|ti] corresponding to each image. Using this data as input, we perform the depth map prediction task and complete the generation of the dense point cloud.

[0028] Step 2.13 is to use the point cloud and normal vector to Poisson reconstruct the three-dimensional surface and obtain its implicit function;

[0029] In step 2.14, after obtaining the implicit function, the Marching Cubes algorithm is used to extract the three-dimensional mesh of the object surface and establish a three-dimensional tunnel model.

[0030] In step 3, specifically:

[0031] Step 3.1, import the established 3D tunnel model into Blender software, and use S to scale and R to rotate to adjust the posture;

[0032] Step 3.2: For camera internal parameter settings, input focal length, sensor size, and image resolution. The principal point coordinates are located at the image geometric center by default.

[0033] Step 3.3, based on the input data of step 3.2, obtain the camera's intrinsic parameter matrix K;

[0034] Step 3.4: For the camera external parameter setting, the camera position is set, and the camera translation vector is determined based on its spatial position and orientation information in the world coordinate system. The camera's specific position translation vector T in three-dimensional space is determined by setting the camera's three-dimensional position (Tx, Ty, Tz).

[0035] In step 3.5, the camera's orientation is achieved by setting its Euler angles (θx, θy, θz), where θx controls the camera's rotation around the X axis; θy controls the camera's rotation around the Y axis; and θz controls the camera's rotation around the Z axis.

[0036] Step 3.6, the Euler angle rotation order in Blender is XYZ Euler, and the rotation matrix is: R = Rx(θx)·Ry(θy)·Rz(θz)

[0037] The three components are the rotation matrices around the X, Y, and Z axes, respectively. The formulas are as follows:

[0038]

[0039] Step 3.7, based on the camera translation vector T in step 3.4 and the camera rotation matrix in step 3.6, the camera posture matrix is ​​obtained as [R|T];

[0040] Step 3.8: For the long linear structure of the tunnel wall, virtual cameras are placed at equal intervals along the tunnel axis; each camera is set with a fixed elevation angle;

[0041] Step 3.9, set the image resolution; when rendering the camera view, export the depth map D(u,v) corresponding to each image;

[0042] Step 3.10, activate each camera object in sequence according to the number and render the current view image;

[0043] Step 3.11: After rendering, the image file is automatically named according to the camera sequence number and saved. The image file name maintains a one-to-one correspondence with the camera layout order.

[0044] Step 3.12, traverse the mesh structure of the current model object, read the local coordinates v.co of all vertices, and use the object's world coordinate transformation matrix to uniformly map them to the global 3D coordinate system. Each vertex coordinate undergoes matrix multiplication; the converted world coordinate points are output as a plain text .txt file, with one 3D point per line.

[0045] Step 3.13: After completing the 3D model layout and virtual camera configuration, derive the 3D world coordinates of the vertices. Based on the known camera intrinsic parameter matrix K and extrinsic parameter matrix [R|T], and the projection depth D(u,v) in the camera coordinate system, project all 3D points onto the image pixel plane, establishing a one-to-one mapping relationship between image pixel coordinates and 3D points.

[0046] In step 3.14, according to the mapping relationship in step 3.13, the three-dimensional point (Xw, Yw, Zw) is mapped to the image pixel coordinates (u, v). By discretizing the pixel coordinates and constructing a hash table structure with (u, v) as the key and the three-dimensional point as the value, a one-to-one mapping between the image and the three-dimensional space point is achieved.

[0047] In step 4, specifically:

[0048] Step 4.1, grayscale processing is performed on the multiple tunnel wall images obtained in step 3.11;

[0049] Step 4.2, apply the Canny edge detection algorithm to the grayscale image to extract the edge structure information in the image;

[0050] Step 4.3: Count the ratio of the number of edge pixels to the total number of pixels in the image as a preliminary indicator of crack presence. When this ratio exceeds a set threshold, the image is preliminarily judged to have crack structures, and the process proceeds to the next step. If the ratio is below the threshold, the image structure is considered intact, with no obvious cracks. Only the image number and the "no crack" status label are recorded, and the process ends.

[0051] Step 4.4: Denoise the cracked image to reduce image noise. Apply bilateral filtering to remove noise while preserving edge details and perform crack enhancement.

[0052] Step 4.5: Apply an adaptive threshold segmentation algorithm to the filtered and enhanced image to automatically separate the crack area from the background;

[0053] Step 4.6: Use the binary image output from step 4.5.5 as input, with the pixel value of the crack area set to 0, i.e., black; the pixel value of the background area set to 255, i.e., white; and perform connected domain analysis on the black pixels in the image using the 8-neighborhood definition method.

[0054] Step 4.7: Calculate the pixel area of ​​each connected domain and set an area threshold. If the area of ​​the connected domain is less than the threshold, it is considered to be unstructured noise and deleted. Otherwise, it is retained and proceeds to the next step.

[0055] Step 4.8: For the connected domain that remains, extract its minimum bounding rectangle and calculate its aspect ratio. Set a threshold. If the aspect ratio of a region is lower than this threshold, it is considered a noise region and removed.

[0056] Step 4.9: For connected domains that satisfy both the area greater than the threshold and the morphological characteristics showing obvious directionality, retain their structure in the binary image;

[0057] Step 4.10: Use Zhang-Suen thinning algorithm to perform skeleton extraction on the crack area;

[0058] Step 4.11: After skeleton extraction, the skeleton points are classified based on the pixel neighborhood topology to identify four types of pixels: isolated points, endpoints, crack points, and intersection points. All isolated points are directly removed; intersection points are used as markers for complex structural areas.

[0059] Step 4.12: Starting from the identified crack endpoints, pixel tracking is performed based on the skeleton connectivity relationship. Complete line segments are extracted and their pixel coordinates are recorded. Their lengths are calculated based on the Euclidean distance. A length threshold is set. If the line segment length is less than the threshold, it is considered an invalid burr and the entire segment is removed.

[0060] Step 4.13: After the two-dimensional pixel point of the crack is determined, the three-dimensional point is calculated by back-projection based on the camera intrinsic parameter matrix K in step 3.3, the extrinsic parameter matrix [R|T] in step 3.7, and the projection depth D(u,v) in step 3.9;

[0061] Step 4.14, based on the mapped three-dimensional coordinates, use the three-dimensional Euclidean distance to quantitatively measure the entire crack trace;

[0062] Step 4.15: Fit the plane equation to each fracture line segment, unify the normal vector direction, sort the planes by the normal vector direction, and calculate the vertical distances between adjacent planes in sequence;

[0063] Step 4.16: After completing the 3D coordinate recovery and crack parameter quantification, based on the image number and pixel coordinates corresponding to each crack, combined with the internal and external parameter matrix and camera posture, the 3D space coordinates are regressed to the overall 3D model for annotation.

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

[0065] The method described in this paper generates a 3D model and uses a virtual camera to render virtual images from multiple angles and positions. Compared to a single 2D image, it can provide higher-precision positioning and structural identification in complex tunnel environments. By combining texture mapping technology with the 3D model, an ideal virtual image is generated, overcoming the perspective, lighting, and spatial limitations of actual image acquisition. This allows for high-precision identification of rock mass structural elements on tunnel walls, providing reliable data support for tunnel design, construction, and subsequent maintenance. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 It is a sparse 3D point cloud obtained by triangulation and pose estimation methods;

[0067] Figure 2 It is the front view of the three-dimensional tunnel modeling graphics;

[0068] Figure 3 It is a side view of the three-dimensional tunnel modeling graphic;

[0069] Figure 4 This is the result of separating the crack area and the background through the adaptive threshold segmentation algorithm;

[0070] Figure 5 This is the result image after applying the connected domain denoising process;

[0071] Figure 6 is the crack skeleton map refined using the Zhang-Suen refinement algorithm;

[0072] Figure 7 This is a diagram of the crack skeleton structure optimized by burr removal. DETAILED DESCRIPTION

[0073] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0074] Example 1

[0075] The present invention is based on a tunnel fissure identification method based on 3D real scene reconstruction and orderly re-collection by a virtual camera, and is specifically implemented according to the following steps:

[0076] Step 1, obtaining rock mass crack images;

[0077] Shoot a video of rock fractures and export the video as a sequence of frames in PNG format, numbering them in sequence for subsequent processing.

[0078] Step 2: 3D modeling is used to obtain 3D point cloud data of the tunnel wall from the collected rock mass fracture images, and the geometry of the tunnel is reconstructed. The acquired texture information is mapped to the 3D mesh vertices using an image processing algorithm to generate a complete 3D model and texture map.

[0079] Specifically:

[0080] Step 2.1: For each rock fracture image, a set of key points with scale and rotation invariance is extracted using the scale-invariant feature transform (SIFT) method, and a corresponding description vector (128 dimensions) is generated for each key point, which is the feature point descriptor used to characterize its local image structure characteristics;

[0081] Step 2.2: For the descriptor sets N×128 and M×128 in the rock mass fracture images A and B of adjacent frames, the feature point matching relationship between the images is calculated based on the similarity between the description vectors, and a 2D-2D corresponding point set between the images is established;

[0082] Specifically, let x be the descriptor of the feature point of rock fracture image A, traverse all descriptors in rock fracture image B, calculate its Euclidean distance with x, and select the two descriptors y1 and y2 of the nearest neighbor and the second nearest neighbor. If the ratio constraint d(x,y1) / d(x,y2)<0.75 is satisfied, the matching is considered valid, and the descriptor y1 corresponding to x is obtained. Finally, the matching point pair set {(p1 (i) ,p2 (i) )};

[0083] The Euclidean distance calculation formula is as follows:

[0084]

[0085] Among them, x i ,y i They are a matching pair to each other;

[0086] Step 2.3: Use the RANSAC algorithm to estimate the essential matrix; specifically:

[0087] According to the pixel coordinates of the rock mass crack images in adjacent frames, the matching point pair set {(p1 (i) ,p2 (i))} and the camera intrinsic parameter matrix K, which converts the image coordinates from pixel coordinates to the normalized coordinate system (camera coordinate system), x1 = K -1 ·p 1 ,x2=K -1 ·p 2 ; x1, x2 are the coordinates of the matching point on the normalized plane (i.e., the plane with Z=1 in the camera coordinate system);

[0088] For each set of matching points (x1 (i) ,x2 (i) ), E is the essential matrix, E∈R 3×3 , then the essential matrix satisfies

[0089] Use the eight-point method (based on linear least squares method) to solve the essential matrix E, that is, select eight pairs of matching feature points in the matching pair set and construct a matrix A such that

[0090] That is, each pair of points constructs a row of A i = = [x 2x x 1x ,x 2x x 1y ,x 2x ,x 2y x 1x ,x 2y x 1y ,x 2y ,x 1x ,x 1y ,1], all A i Composition matrix Where e is the column vector of the essential matrix E;

[0091] Perform singular value decomposition (SVD) on the matrix A, take the singular vector corresponding to the minimum singular value as the solution of e, and reconstruct it into a 3×3 matrix, which is the preliminary estimated essential matrix E1;

[0092] Step 2.4, perform singular value decomposition on the initially estimated essential matrix E1 to obtain the estimated essential matrix E2 = UDV T , where U and V are orthogonal matrices and D is a diagonal matrix. Then, to satisfy the theoretical requirement that the essential matrix rank is 2, the third singular value in D is set to zero, and the modified matrix D′=diag(s,s,0) is constructed and the essential matrix E′=UD′V is reconstructed. T At the same time, the U and V matrices are retained for posture recovery in subsequent steps;

[0093] Step 2.5, according to the formula Calculate the epipolar geometric error of all point pairs. When the error is less than the threshold of 1.5, the point is marked as an "inlier point", otherwise it is an "outlier point". After several rounds of iteration, the E2 with the most inliers is retained as the optimal solution. Output the estimated essential matrix E2 and the corresponding inlier point set.

[0094] Step 2.6: Based on the modified essential matrix E′=UD′V obtained in step 2.4 T , use the U and V matrices to decompose the relative posture of the camera (rotation matrix R and translation vector t), and construct four sets of candidate (R, t) combinations according to the standard decomposition method. Each set of rotation and translation corresponds to a different camera posture. The rotation selection matrix Rz represents the standard orthogonal matrix of ±π / 2 rotation around the z axis;

[0095] R1=UR Z (π / 2)V T ,t1=U[:,2]

[0096] R2=UR Z (-π / 2)V T ,t2=U[:,2]

[0097] R3=UR Z (π / 2)V T ,t3=-U[:,2]

[0098] R4=UR Z (-π / 2)V T ,t4=-U[:,2]

[0099] In step 2.7, only one of the four solutions is correct. That is, under this combination, all the 3D points obtained by triangulation of the matching points are simultaneously in front of both cameras (i.e., the Z coordinate is positive). The correct solution is determined by checking the 3D points, thus obtaining the pose matrix P = [R|t] for the current frame.

[0100] Step 2.8: Use triangulation to calculate the coordinates of the 3D points. First, calculate the projection matrix M for each rock fracture image based on the known intrinsic parameter matrix K and extrinsic parameter matrix (rotation matrix R and translation vector t) of the camera. Let u be the homogeneous coordinates of the 2D point, U be the homogeneous coordinates of the 3D point, and M be the projection matrix. Then the projection matrix satisfies u = MU. The projection matrix can be solved using the formula:

[0101] M=K(R|t)

[0102] Step 2.9, based on the projection matrix and corresponding points of the rock fracture images of adjacent frames, the projection matrix of the i-th camera is M i =K i [R i ,t i ]=[M i1 ,Mi2 ,M i3 ] T , the corresponding pixel coordinate is x i =[u i ,v i ,1] T According to the projection equation, we can get d i x i =M i X, cross-multiply both sides of the equation by x i You can get x i ×M i X = 0. Rearranging the equations yields the scene's three-dimensional point cloud X;

[0103]

[0104] In step 2.10, the four candidate poses are verified according to the positive depth constraint (i.e., the Z coordinates of the three-dimensional points obtained by triangulation are all positive in the two camera coordinate systems). Finally, the one that satisfies the positive Z value of all points is selected as the camera pose [R|t] of the current frame image, which serves as the basis for subsequent triangulation and error optimization.

[0105] Step 2.11, after obtaining preliminary 3D points through triangulation, perform bundle adjustment optimization to minimize the reprojection error of all 3D points and camera parameters. The error function is

[0106]

[0107] Among them, u ij is the pixel coordinate, M j is the camera projection matrix, X i are the three-dimensional point coordinates;

[0108] During the optimization process, the LM (Levenberg-Marquardt) nonlinear least squares algorithm is used for iterative solution. By minimizing the reprojection error, the optimized updated 3D point cloud coordinates and the updated pose matrix image corresponding to each frame image are obtained.

[0109] Step 2.12: Based on the acquired sparse 3D point cloud X, the camera pose matrix [R|t], and the intrinsic parameter matrix K, we further calculate the projection matrix Mi = K[Ri|ti] corresponding to each image. Using this data as input, we perform the depth map prediction task and complete the generation of the dense point cloud.

[0110] Specifically: First, for each image, according to different depth values, set d∈{d1,d2,...,d N}, use the homography transformation formula to transform the features of other images into the coordinate system of this image, calculate the pixel coordinates u′ on the reference image after projection, and construct the cost volume to represent the degree of matching of each pixel in depth. The formula is:

[0111] u′=Hx=K′[R-tn T / d]K -1 u

[0112] Where H is the homography matrix, K and K′ are the intrinsic parameter matrices of the two images, R and t are the rotation and translation matrices between two adjacent frames, n is the normal vector of the plane scan, and d is the depth value of the plane scan;

[0113] After completing the projection of feature points of all adjacent images and completing the preliminary alignment, according to the set multiple depth hypotheses, each reference image pixel (x, y) is k Based on the transformed projection position u′, the corresponding pixel values ​​are sampled in multiple perspective images, and the difference between them and the pixel values ​​in the current reference image is calculated, and this is used as the matching cost to finally construct the three-dimensional cost volume C(x, y, d k ), which is used to describe the matching consistency of the pixel at different depths.

[0114] The cost body construction formula is:

[0115]

[0116] Where: I0 is the current reference image, I j is the jth source image; u′, v′ are the pixels according to the depth d k and the corresponding pixel coordinates after projection using the homography transformation formula; W represents the local window size (5×5), which is used to aggregate pixel differences within the matching area; M is the number of source images involved in constructing the cost.

[0117] For each pixel (x, y) in the reference image, based on its corresponding candidate 3D point (triangulated by multiple depth d hypotheses), calculate the triangulation, resolution, and angle of incidence between it and each auxiliary image. On the constructed cost volume C(x, y, d), the depth with the minimum cost is selected as the initial estimate:

[0118] Calculate the depth value that best meets the photometric consistency assumption.

[0119] The estimated depth value d *(x, y), back-project the pixel into three-dimensional space to obtain its three-dimensional coordinate point X(x, y). Select the 5×5 local neighborhood of the pixel in the image and obtain the back-projection point set within its neighborhood. Construct the local normal vector n(x, y) by performing a least squares plane fitting operation on these neighborhood points. Use principal component analysis (PCA) and take the eigenvector corresponding to its minimum eigenvalue as the normal direction. Finally, output the unit normal vector n(x, y) for each pixel position;

[0120] Finally, the PatchMatch algorithm is used to iteratively update the depth and normal vectors until the convergence conditions are met or the preset number of iterations is reached, and a high-precision dense three-dimensional point cloud is output.

[0121] Step 2.13, Poisson reconstruction uses point clouds and normal vectors to reconstruct a 3D surface. There exists an implicit function whose gradient field is equal to or close to the normal field of the point cloud. This implicit function can be considered an indicator function, which is 1 inside the object, 0 outside the object, and 0.5 on the surface of the object. The general form of the Poisson equation is:

[0122]

[0123] Where Δ is the Laplace operator, is the gradient operator, χ is the implicit function to be solved, is the normal vector field of the point cloud;

[0124] Step 2.14: After obtaining the implicit function, use the Marching Cubes algorithm to extract the three-dimensional mesh of the object surface and build a three-dimensional tunnel model;

[0125] Step 3: Use Blender software to set the parameters of the virtual camera, obtain the captured virtual image, and match the position of the virtual camera with the real position of the model;

[0126] Step 3.1, import the 3D tunnel model (.obj format) into Blender. After importing, press G to move the model to the appropriate position, and use S to scale and R to rotate to adjust the posture.

[0127] In step 3.2, set the camera's internal parameters by inputting the focal length, sensor size, and image resolution. The principal point coordinates are located at the image's geometric center by default.

[0128] Step 3.3, based on the input data of step 3.2, obtain the camera's intrinsic parameter matrix K;

[0129]

[0130] in,

[0131] fx, fy are the focal lengths (horizontal and vertical directions) at the pixel scale;

[0132] u0, v0 are the main points (image center coordinates);

[0133] f is the focal length, in mm;

[0134] image_width is the image width, in px (pixels);

[0135] image_height is the image height, in px (pixels);

[0136] sensor_width is the sensor width in mm;

[0137] Step 3.4, for the camera external parameter setting, depends on its spatial position and orientation information in the world coordinate system, set the camera position, determine the camera translation vector, and determine the camera's specific position translation vector T in three-dimensional space by setting the camera's three-dimensional position (Tx, Ty, Tz):

[0138] T=[T X T Y T Z ] T

[0139] Among them, Tx, Ty, Tz are the camera translation amounts, indicating the camera's position in the world coordinate system;

[0140] In step 3.5, the camera's orientation is achieved by setting its Euler angles (θx, θy, θz), where:

[0141] θx is the rotation angle of the camera around the X axis (pitch angle), which is used to ensure that the lens is aligned with the center of the structure or wall;

[0142] θy controls the rotation around the Y axis (yaw angle) and adjusts the left and right viewing angles of the camera;

[0143] θz controls the rotation around the Z axis (roll angle) to ensure the image is horizontal.

[0144] Step 3.6, the Euler angle rotation order in Blender is XYZ Euler, and the rotation matrix is:

[0145] R=Rx(θx)·Ry(θy)·Rz(θz)

[0146] The three components are the rotation matrices around the X, Y, and Z axes, respectively. The formulas are as follows:

[0147]

[0148] Step 3.7, based on the camera translation vector T in step 3.4 and the camera rotation matrix in step 3.6, the camera posture matrix is ​​obtained as [R|T];

[0149] Step 3.8: For the long linear structure of the tunnel wall, in Blender, use Python script to batch-space virtual cameras at equal intervals along the tunnel axis. Each camera is set with a fixed pitch angle (e.g., rotated 90° around the X axis) to ensure that the lens faces the tunnel wall. The position of the virtual camera is defined as follows: the spatial position of the i-th virtual camera is defined as T i =[T x ,T y ,T z +i·Δz]T,i=0,1,...,N;

[0150] Where Δz is the layout spacing, and N is the total number of deployed cameras.

[0151] In step 3.9, set the image resolution and output image format to PNG; select Eevee or Cycles as the rendering engine, set the ambient lighting (World Strength recommended value 1-2) and enable the exposure control function (exposure value recommended value +1.0), and turn off the motion blur effect. When rendering the camera view, select "Z" output and export the depth map D(u,v) corresponding to each image.

[0152] Step 3.10: Call the Blender rendering interface through the Python script, activate each camera object in sequence according to the number, set it as the current camera, and call the command bpy.ops.render.render(write_still=True) to render the current perspective image.

[0153] Step 3.11: After rendering, the image file is automatically named according to the camera serial number. The save path can be preset by the user. The image file name maintains a one-to-one correspondence with the camera layout order.

[0154] In step 3.12, use Blender's built-in bpy scripting interface to traverse the mesh structure of the current model object (Object), read the local coordinates v.co of all vertices, and use the object's world coordinate transformation matrix obj.matrix_world to uniformly map them to the global 3D coordinate system. Each vertex coordinate undergoes matrix multiplication:

[0155] X world =M world ·X local

[0156] Among them, Mworld is the world matrix of the object obj.matrix_world;

[0157] Xlocal is the model vertex v.co.

[0158] The converted world coordinate points are output in the form of a plain text .txt file, with each line recording a three-dimensional point (Xw, Yw, Zw), providing input data for establishing the image space-three-dimensional space mapping relationship.

[0159] In step 3.13, after completing the 3D model layout and virtual camera configuration in Blender, export the 3D world coordinates of the vertices. Using a Python script, we project all 3D points onto the image pixel plane based on the known camera intrinsic and extrinsic matrix [R|T], as well as the projection depth D(u,v) in the camera coordinate system. This establishes a one-to-one mapping between image pixel coordinates and 3D points. The specific mapping relationship is as follows:

[0160]

[0161] In step 3.14, following the mapping relationship from step 3.13, map the 3D point (Xw, Yw, Zw) to the image pixel coordinates (u, v). This is done by discretizing the pixel coordinates and constructing a hash table structure with (u, v) as the key and the 3D point as the value, achieving a one-to-one mapping between the image and the 3D point. This structure is stored in JSON format, supporting fast forward index queries.

[0162] Step 4: Use geometric methods to extract the geometric information of the crack (such as length, width, direction, etc.) from the acquired virtual image, and project the two-dimensional image coordinates back into the three-dimensional space through the projection matrix to obtain the spatial coordinates of the crack.

[0163] Example 2

[0164] Further, step 4 is specifically as follows:

[0165] Step 4.1: Preprocess the multiple tunnel wall images obtained in step 3.11. First, grayscale the images to eliminate color interference and unify the image channel format. After loading the original color image, convert it to a grayscale image using the weighted average formula:

[0166] Im=0.299·I R +0.587·I G +0.114·I B

[0167] Among them, I R ,I G ,I B are the pixel values ​​of red, green and blue channels respectively, and Im is the grayscale image;

[0168] Step 4.2: Apply the Canny edge detection algorithm to the grayscale image to extract the edge structure information in the image in order to identify the possible crack edge contours;

[0169] Step 4.3: Count the ratio of the number of edge pixels to the total number of pixels in the image, R = N edge \N total , as a preliminary indicator of crack existence. When the ratio R exceeds the set threshold (e.g., 2%), it is preliminarily determined that there is a crack structure in the image, and the next step is entered; if the ratio is lower than the threshold, the image structure is considered intact, with no obvious cracks, and only the image number and the "no crack" status label are recorded, and the process ends;

[0170] In step 4.4, denoising is performed on the crack image to reduce image noise, and bilateral filtering is applied to remove noise while retaining edge details to enhance cracks.

[0171] Bilateral filtering formula:

[0172]

[0173] Among them, G d is the spatial adjacency Gaussian function, G r is the Gaussian function of pixel value similarity;

[0174] i: represents the index of the center pixel position currently being processed;

[0175] j: represents the position index of other pixels in the neighborhood of pixel i;

[0176] I i , I j Respectively represent the grayscale values ​​of pixels at positions i and j in the image;

[0177] N represents the set of neighboring pixels of pixel i;

[0178] ω i : Normalization factor, representing the sum of all weights, used to maintain the consistency of brightness after filtering;

[0179] δ BF (I i ): represents the new value of the pixel at position i after bilateral filtering.

[0180] Step 4.5: Apply an adaptive threshold segmentation algorithm to the filtered and enhanced image to automatically separate the crack area from the background. Specifically:

[0181] Step 4.5.1: Select the local neighborhood area of ​​any pixel (x, y) in the image for feature analysis and define a k×k square sampling window. The window size is usually set to an odd number.

[0182] Step 4.5.2: Construct a Gaussian window, where (i, j) is the point coordinate, and the center position is used as the origin for sampling. The Gaussian weight is calculated as follows:

[0183]

[0184] In step 4.5.3, for each pixel (x, y), calculate the local Gaussian weighted average according to the following formula:

[0185]

[0186] Among them, I(x+i,y+j) is the pixel value in the neighborhood window, and w(i,j) is the corresponding Gaussian weight.

[0187] In step 4.5.4, calculate the local adaptive threshold for each pixel using the following formula:

[0188] T(x,y)=mean(x,y)-C

[0189] Where C is a constant subtracted from the Gaussian weighted mean, which is usually used to adjust the sensitivity of the threshold.

[0190] Step 4.5.5: Perform image threshold segmentation according to the local threshold segmentation rule:

[0191]

[0192] Compare the size of the pixel point with the local adaptive threshold. If it is less than the threshold, it is segmented as a crack target (black); if it is greater than or equal to the threshold, it is segmented as the background (white).

[0193] In step 4.6, the binary image output from step 4.5.5 is used as input. The pixel value of the crack area is 0 (black), and the pixel value of the background area is 255 (white). The connected domain analysis is performed on the black pixels (value 0) in the image using the 8-neighborhood definition method.

[0194] Step 4.7, calculate the pixel area of ​​each connected domain (i.e. the number of black pixels), set the area threshold (such as A min =35 pixels). If the area of ​​a connected domain is smaller than the threshold, it is considered as unstructured noise and is deleted (the area is set to background pixels 255); otherwise, it is retained and proceeds to the next step of judgment;

[0195] Step 4.8: For the connected domain that remains, extract its minimum bounding rectangle and calculate its aspect ratio γ = L\W, where L is the longer side and W is the shorter side. Set the discrimination threshold (e.g. γ min =1.2), if the aspect ratio of a region is lower than this value (i.e. the shape is approximately circular), it is considered as a noise region and is removed;

[0196] In step 4.9, for connected domains that satisfy both the area greater than the threshold and the morphological features showing obvious directionality, their structures in the binary image are retained and used as the input image for subsequent skeleton extraction and crack parameter quantification;

[0197] In step 4.10, the Zhang-Suen thinning algorithm is used to perform skeleton extraction on the crack area. This algorithm is based on an 8-neighborhood and sequentially scans all black pixels in the image, performing pixel culling in two iterative steps:

[0198] The first step of elimination conditions includes:

[0199] The number of black pixels in the neighborhood is 2≤N(P0)≤6;

[0200] Number of changes from white to black A(P1) = 1;

[0201] The structural conditions P1·P3·P5=0 and P3·P5·P7=0 are maintained;

[0202] The second exclusion criteria include:

[0203] The number of black pixels in the neighborhood is 2≤N(P0)≤6;

[0204] Number of changes from white to black A(P1) = 1;

[0205] The structural conditions P1·P5·P7=0 and P1·P3·P7=0 are maintained;

[0206] The algorithm gradually removes boundary pixels through iterative operations until there are no more deletable pixels in the image.

[0207] In step 4.11, after skeleton extraction, some cracks may still have endpoint glitches, isolated pixels, or redundant short branches. To improve the structural purity, the skeleton points are classified based on the pixel neighborhood topology and four types of pixels are identified:

[0208] Isolated point: the number of black pixels in the neighborhood n = 0;

[0209] Endpoint: number of black pixels in the neighborhood n = 1;

[0210] Crack points: n = 2 and the adjacency relationship is linearly continuous;

[0211] Crossover point: n ≥ 3;

[0212] All isolated points will be directly removed; intersection points are used as markers for complex structural areas, and the endpoints record coordinates for subsequent structure tracking.

[0213] Step 4.12: Starting from the identified crack endpoint, perform pixel tracking according to the skeleton connectivity relationship, extract the complete line segment and record its pixel coordinate set {a1, a2, ..., a m}, calculate its length L according to the Euclidean distance;

[0214]

[0215] Among them, x k ,y k Represents the horizontal and vertical coordinates of the k-th pixel point in the skeleton segment in the image coordinate system;

[0216] Set the length threshold (such as L min =9 pixels), if the length of a line segment is less than the threshold, it is determined to be an invalid burr and the entire segment is removed.

[0217] Step 4.13: After preprocessing the image and extracting the cracks, the 2D pixel points of the cracks are determined. Based on the camera intrinsic parameter matrix K from step 3.3, the extrinsic parameter matrix [R|T] from step 3.7, and the projection depth D(u,v) from step 3.9, the 3D point (Xw, Yw, Zw) is calculated by back-projection:

[0218]

[0219] In step 4.14, based on the mapped 3D coordinates, the entire crack trace is quantitatively measured using the 3D Euclidean distance. The crack length calculation formula is as follows:

[0220]

[0221] Where l represents the actual extension length of the crack trace, n corresponds to the number of three-dimensional coordinate points contained in a single trace, (x i ,y i ,z i ) represents the three-dimensional coordinates of the i-th point.

[0222] Step 4.15: Fit each fracture line segment with the plane equation Ax+By+Cz+D=0. After unifying the normal vector direction, sort the planes by the normal vector direction and calculate the vertical distances between adjacent planes in sequence.

[0223]

[0224] In step 4.16, after completing 3D coordinate recovery and crack parameter quantification, based on the image number and pixel coordinates (u, v) corresponding to each crack, combined with the internal and external parameter matrices and camera pose, the 3D space coordinates (Xw, Yw, Zw) are regressed to the overall 3D model for annotation in the Blender Python script. The model is then exported to .obj format and displayed superimposed on the original model.

[0225] Example 3

[0226] Figure 1 This is a sparse 3D point cloud generated through the triangulation and pose estimation methods in steps 2.13 to 2.17. The 3D points in the image represent the calculated 3D coordinates of matching points obtained from different camera viewpoints. This point cloud image demonstrates the feasibility of recovering a 3D scene from 2D image information and provides reliable foundational data for subsequent dense point cloud generation and 3D reconstruction.

[0227] Figure 2 and Figure 3 These are the front and side views of the 3D tunnel modeling graphics verified in step 2. Key information such as the texture, cracks, and structural features of the rock surface can be observed. The precise modeling method provides a basis for the subsequent provision of crack images.

[0228] Example 4

[0229] Table 1, Step 3: From importing the 3D tunnel model and placing the virtual camera to image rendering and coordinate mapping, the corresponding pixel coordinates are obtained. Ultimately, the mapping relationship between these pixel coordinates and world coordinates is stored in a hash table for management. This hash table structure not only ensures efficient data storage but also verifies that the correspondence between the image and the actual world coordinates is reliable and traceable, from virtual camera rendering to 3D reconstruction.

[0230] Table 1 Hash table storage results

[0231] Hash table index Key-value pair linked list Hash table index Key-value pair linked list index=270007 (270,7):(2.83,3.17,3.82) index=317232 (317,232):(2.76,3.24,3.35) index=290031 (290,31):(2.80,3.20,3.77) index=324290 (324,290):(2.75,3.25,3.23) index=297055 (297,55):(2.79,3.21,3.72) index=327314 (327,314):(2.74,3.25,3.18) index=303088 (303,88):(2.78,3.22,3.65) index=330333 (330,333):(2.74,3.26,3.14) index=303131 (303,131):(2.78,3.22,3.56) index=364395 (364,395):(2.69,3.31,3.01) index=310184 (310,184):(2.77,3.23,3.45) index=419481 (419,481):(2.61,3.39,2.83)

[0232] Example 5

[0233] Figure 4-Figure 7 This is the change process of crack image processing in step 4. In the crack extraction process, the crack area and background are separated by the adaptive threshold segmentation algorithm in step 4.5 to accurately extract the crack information ( Figure 4 Then, steps 4.6-4.9 apply connected domain denoising to remove noise and retain the main crack structure ( Figure 5 ). Subsequently, step 4.10 uses the Zhang-Suen refinement algorithm to refine the crack skeleton and extract a clear centerline ( Figure 6Finally, step 4.11 optimizes the crack skeleton structure by removing burrs, removing redundant points and short branches ( Figure 7 The entire process allows the cracks to be gradually extracted from the background and refined, ensuring accurate extraction and optimization of the cracks.

[0234] Example 6

[0235] The tunnel fissure identification method based on 3D real scene reconstruction and orderly re-collection by virtual camera of the present invention has the following advantages:

[0236] (1) The cracks in the virtual image are automatically identified through image processing technology, and the cracks are projected back from the two-dimensional image to the three-dimensional model in combination with the projection matrix to restore the spatial coordinates of the structural features and ensure the accurate positioning and quantitative measurement of the rock structure elements in three-dimensional space.

[0237] (2) Through three-dimensional modeling, images from multiple perspectives are spliced ​​into a complete and continuous tunnel wall model, effectively restoring the actual spatial structure and avoiding problems such as structural fracture and recognition fragmentation caused by occlusion or incomplete perspective in traditional image recognition. This enables continuous recognition and complete extraction of cracks, significantly improving recognition integrity.

[0238] The present invention's tunnel fracture identification method, based on 3D real-scene reconstruction and sequential re-acquisition with a virtual camera, utilizes the SFM and MVS algorithms to generate a 3D model with texture mapping. It then employs virtual camera technology to render high-quality images from multiple angles within the Blender environment. Image processing and feature extraction algorithms automatically identify fracture parameters. Using a projection matrix, the 2D image coordinates are back-projected into 3D space to restore their spatial position. This method quantitatively extracts fracture structural elements and labels their spatial positions, achieving high-precision identification of complex rock mass fractures and providing crucial data support for tunnel stability analysis and structural modeling.

Claims

1. A tunnel fissure identification method based on 3D real scene reconstruction and virtual camera sequential re-collection is characterized by: Please follow the steps below to implement it: Step 1, obtaining rock mass crack images; Step 2: Based on the collected rock fracture images, obtain 3D point cloud data of the tunnel wall, reconstruct the tunnel geometry, and map the acquired texture information to 3D mesh vertices using an image processing algorithm to generate a 3D tunnel model. Step 3: Set the parameters of the virtual camera to obtain the captured virtual image and match the position of the virtual camera with the real position of the model; Step 4: Extract the information of the cracks in the virtual image and project the two-dimensional image coordinates back into the three-dimensional space through the projection matrix to obtain the spatial coordinates of the cracks.

2. The method for identifying tunnel fissures based on 3D real scene reconstruction and orderly re-acquisition of virtual cameras according to claim 1, characterized in that: In the step 1, specifically: shoot a video of rock fractures, and export the video as a sequence of frames in PNG format, and number them in sequence.

3. The tunnel fissure identification method based on 3D real scene reconstruction and virtual camera orderly re-acquisition according to claim 1, characterized in that: In the step 2, specifically: Step 2.1: For each rock fracture image, use the SIFT method to extract a set of key points with scale and rotation invariance, and generate a corresponding description vector for each key point, which is the feature point descriptor; Step 2.2: For the descriptor sets in the rock mass fracture images A and B in adjacent frames, the feature point matching relationship between the images is calculated based on the similarity between the description vectors, and a 2D-2D corresponding point set between the images is established; Step 2.3, use the RANSAC algorithm to estimate the essential matrix E1; Step 2.4, perform singular value decomposition on the estimated essential matrix E1 to obtain the estimated essential matrix E2 = UDV T , where U and V are orthogonal matrices, D is a diagonal matrix, and the third singular value in D is set to zero, constructing the modified matrix D′=diag(s,s,0) and reconstructing the essential matrix E′=UD′V T At the same time, keep the U and V matrices; Step 2.5, according to the formula Calculate the epipolar geometric error of all point pairs. If the error is less than the threshold of 1.5, the point is marked as an "inlier point", otherwise it is an "outlier point". After several rounds of iteration, the E2 with the most inliers is retained as the optimal solution. Output the estimated essential matrix E2 and the corresponding inlier point set. Step 2.6: Based on the modified essential matrix E′=UD′V obtained in step 2.4 T , use the U and V matrices to decompose the relative posture of the camera, and construct four sets of candidate (R, t) combinations according to the standard decomposition method. Each set of rotation and translation corresponds to a different camera posture. The rotation selection matrix Rz represents the standard orthogonal matrix of ±π / 2 rotation around the z axis; R1=UR Z (π / 2)V T ,t1=U[:,2] R2=UR Z (-π / 2)V T ,t2=U[:,2] R3=UR Z (π / 2)V T ,t3=-U[:,2] R4=UR Z (-π / 2)V T ,t4=-U[:,2] In step 2.7, only one of the four solutions is correct, that is, under this combination, all the 3D points obtained by triangulation of the matching points are simultaneously in front of both cameras. The correct solution is determined by checking the 3D points, thus obtaining the pose matrix P = [R|t] for the current frame. Step 2.8, use triangulation to calculate the coordinates of the three-dimensional points. First, calculate the projection matrix M of each rock fracture image based on the known intrinsic parameter matrix K and extrinsic parameter matrix of the camera; Let u be the homogeneous coordinates of a two-dimensional point, U be the homogeneous coordinates of a three-dimensional point, and M be the projection matrix. Then the projection matrix satisfies u=MU; Step 2.9, based on the projection matrix and corresponding points of the rock fracture images of adjacent frames, the projection matrix of the i-th camera is M i =K i [R i ,t i ]=[M i1 ,M i2 ,M i3 ] T , the corresponding pixel coordinate is x i =[u i ,v i ,1] T According to the projection equation, we can get d i x i =M i X, cross-multiply both sides of the equation by x i You can get x i ×M i X=0; rearranging the equations yields the scene's three-dimensional point cloud X; Step 2.10: Verify the four candidate poses according to the positive depth constraint, and finally select the one that satisfies all points with positive Z values ​​as the camera pose [R|t] of the current frame image. Step 2.11: After obtaining preliminary 3D points through triangulation, perform bundle adjustment optimization to minimize the reprojection error of all 3D points and camera parameters. Step 2.12: Based on the acquired sparse 3D point cloud X, the camera pose matrix [R|t], and the intrinsic parameter matrix K, we further calculate the projection matrix Mi = K[Ri|ti] corresponding to each image. Using this data as input, we perform the depth map prediction task and complete the generation of the dense point cloud. Step 2.13 is to use the point cloud and normal vector to Poisson reconstruct the three-dimensional surface and obtain its implicit function; In step 2.14, after obtaining the implicit function, the Marching Cubes algorithm is used to extract the three-dimensional mesh of the object surface and establish a three-dimensional tunnel model.

4. The method for identifying tunnel fissures based on 3D real scene reconstruction and orderly re-acquisition of virtual cameras according to claim 3, characterized in that: In the step 3, specifically: Step 3.1, import the established 3D tunnel model into Blender software, and use S to scale and R to rotate to adjust the posture; Step 3.2: For camera internal parameter settings, input focal length, sensor size, and image resolution. The principal point coordinates are located at the image geometric center by default. Step 3.3, based on the input data of step 3.2, obtain the camera's intrinsic parameter matrix K; Step 3.4: For the camera external parameter setting, the camera position is set, and the camera translation vector is determined based on its spatial position and orientation information in the world coordinate system. The camera's specific position translation vector T in three-dimensional space is determined by setting the camera's three-dimensional position (Tx, Ty, Tz). In step 3.5, the camera's orientation is achieved by setting its Euler angles (θx, θy, θz), where θx controls the camera's rotation around the X axis; θy controls the camera's rotation around the Y axis; and θz controls the camera's rotation around the Z axis. Step 3.6, the Euler angle rotation order in Blender is XYZ Euler, and the rotation matrix is: R = Rx(θx)·Ry(θy)·Rz(θz) The three components are the rotation matrices around the X, Y, and Z axes, respectively. The formulas are as follows: Step 3.7, based on the camera translation vector T in step 3.4 and the camera rotation matrix in step 3.6, the camera posture matrix is ​​obtained as [R|T]; Step 3.8: For the long linear structure of the tunnel wall, virtual cameras are placed at equal intervals along the tunnel axis; each camera is set with a fixed elevation angle; Step 3.9, set the image resolution; when rendering the camera view, export the depth map D(u,v) corresponding to each image; Step 3.10, activate each camera object in sequence according to the number and render the current view image; Step 3.11: After rendering, the image file is automatically named according to the camera sequence number and saved. The image file name maintains a one-to-one correspondence with the camera layout order. Step 3.12, traverse the mesh structure of the current model object, read the local coordinates v.co of all vertices, and use the object's world coordinate transformation matrix to uniformly map them to the global 3D coordinate system. Each vertex coordinate undergoes matrix multiplication; the converted world coordinate points are output as a plain text .txt file, with one 3D point per line. Step 3.13: After completing the 3D model layout and virtual camera configuration, derive the 3D world coordinates of the vertices. Based on the known camera intrinsic parameter matrix K and extrinsic parameter matrix [R|T], and the projection depth D(u,v) in the camera coordinate system, project all 3D points onto the image pixel plane, establishing a one-to-one mapping relationship between image pixel coordinates and 3D points. In step 3.14, according to the mapping relationship in step 3.13, the three-dimensional point (Xw, Yw, Zw) is mapped to the image pixel coordinates (u, v). By discretizing the pixel coordinates and constructing a hash table structure with (u, v) as the key and the three-dimensional point as the value, a one-to-one mapping between the image and the three-dimensional space point is achieved.

5. The method for identifying tunnel fissures based on 3D real scene reconstruction and sequential re-acquisition of virtual cameras according to claim 4, characterized in that: In the step 4, specifically: Step 4.1, grayscale processing is performed on the multiple tunnel wall images obtained in step 3.11; Step 4.2, apply the Canny edge detection algorithm to the grayscale image to extract the edge structure information in the image; Step 4.3: Count the ratio of the number of edge pixels to the total number of pixels in the image as a preliminary indicator of crack presence. When this ratio exceeds a set threshold, the image is preliminarily judged to have cracks, and the process proceeds to the next step. If the ratio is below the threshold, the image is considered intact, with no obvious cracks. Only the image number and the "no crack" status label are recorded, and the process ends. Step 4.4: Denoise the cracked image to reduce image noise. Apply bilateral filtering to remove noise while preserving edge details and perform crack enhancement. Step 4.5: Apply an adaptive threshold segmentation algorithm to the filtered and enhanced image to automatically separate the crack area from the background; Step 4.6: Use the binary image output from step 4.5.5 as input, with the pixel value of the crack area set to 0, i.e., black; the pixel value of the background area set to 255, i.e., white; and perform connected domain analysis on the black pixels in the image using the 8-neighborhood definition method. Step 4.7: Calculate the pixel area of ​​each connected domain and set an area threshold. If the area of ​​the connected domain is less than the threshold, it is considered to be unstructured noise and deleted. Otherwise, it is retained and proceeds to the next step. Step 4.8: For the connected domain that remains, extract its minimum bounding rectangle and calculate its aspect ratio; Set the discrimination threshold. If the aspect ratio of the region is lower than this value, it will be considered as a noise region and removed. Step 4.9: For connected domains that satisfy both the area greater than the threshold and the morphological characteristics showing obvious directionality, retain their structure in the binary image; Step 4.10: Use Zhang-Suen thinning algorithm to perform skeleton extraction on the crack area; Step 4.11: After skeleton extraction, the skeleton points are classified based on the pixel neighborhood topology to identify four types of pixels: isolated points, endpoints, crack points, and intersection points. All isolated points will be directly eliminated. Intersections serve as markers for structurally complex areas; Step 4.12: Starting from the identified crack endpoint, perform pixel tracking based on the skeleton connectivity relationship, extract the complete line segment and record its pixel coordinate set, and calculate its length based on the Euclidean distance; Set a length threshold. If the line segment length is less than the threshold, it is considered an invalid burr and the entire segment is removed. Step 4.13: After the two-dimensional pixel point of the crack is determined, the three-dimensional point is calculated by back-projection based on the camera intrinsic parameter matrix K in step 3.3, the extrinsic parameter matrix [R|T] in step 3.7, and the projection depth D(u,v) in step 3.

9. Step 4.14, based on the mapped three-dimensional coordinates, use the three-dimensional Euclidean distance to quantitatively measure the entire crack trace; Step 4.15: Fit the plane equation to each fracture line segment, unify the normal vector direction, sort the planes by the normal vector direction, and calculate the vertical distances between adjacent planes in sequence; Step 4.16: After completing the 3D coordinate recovery and crack parameter quantification, based on the image number and pixel coordinates corresponding to each crack, combined with the internal and external parameter matrix and camera posture, the 3D space coordinates are regressed to the overall 3D model for annotation.

6. The method for identifying tunnel fissures based on 3D real scene reconstruction and sequential re-acquisition of virtual cameras according to claim 5, characterized in that: In the step 4.5, specifically: Step 4.5.1: Select the local neighborhood area of ​​any pixel in the image for feature analysis and define a k×k square sampling window. The window size is usually set to an odd number. Step 4.5.2, construct the Gaussian window; Step 4.5.3, for each pixel, calculate the local Gaussian weighted average; Step 4.5.4, calculate the local adaptive threshold for each pixel; Step 4.5.5, perform image threshold segmentation according to the local threshold segmentation rule; Compare the size of the pixel point with the local adaptive threshold. If it is less than the threshold, it is segmented as a crack target; if it is greater than or equal to the threshold, it is segmented as background.

7. The method for identifying tunnel fissures based on 3D real scene reconstruction and sequential re-acquisition of virtual cameras according to claim 6, characterized in that: In the step 4.10, specifically: Based on the 8-neighborhood, all black pixels in the image are scanned in sequence, and pixel culling is performed iteratively in two steps: The first step of elimination conditions includes: The number of black pixels in the neighborhood is 2≤N(P0)≤6; Number of changes from white to black A(P1) = 1; The structural conditions P1·P3·P5=0 and P3·P5·P7=0 are maintained; The second exclusion criteria include: The number of black pixels in the neighborhood is 2≤N(P0)≤6; Number of changes from white to black A(P1) = 1; The structural conditions P1·P5·P7=0 and P1·P3·P7=0 are maintained; The boundary pixels are gradually removed through iterative operations until there are no more deletable pixels in the image.

8. The method for identifying tunnel fissures based on 3D real scene reconstruction and sequential re-acquisition with a virtual camera according to claim 6, characterized in that: In step 4.14, the crack length calculation formula is as follows: Where l represents the actual extension length of the crack trace, n corresponds to the number of three-dimensional coordinate points contained in a single trace, (x i ,y i ,z i ) represents the three-dimensional coordinates of the i-th point.

Images