Lung bronchus voxel model optimization method for monocular slam dynamic modeling

Abstract

The application provides a lung bronchus voxel model optimization method for monocular SLAM dynamic modeling, which comprises the following steps: replacing the cavity center line with the key frame trajectory of a lung interventional surgery robot estimated by monocular SLAM; then extracting a bronchus voxel model radial section perpendicular to the robot trajectory at intervals of the minimum voxel edge length, and establishing a projection model of the model in the radial direction of the bronchus; then searching the occupied grid on the section based on a radial line segment search method, proposing a curvature consistency graphic interpolation method, and realizing the elimination of the redundant voxels and the filling of the surface pores of the bronchus model in the radial direction. Finally, based on the mapping relationship between the two-dimensional pixels and the three-dimensional voxel units, the section data is converted into a voxel map, and the optimized lung bronchus voxel model is obtained. The application can solve the problems of voxel redundancy and surface pores in the voxel model.

Description

TECHNICAL FIELD

[0001] The application relates to the technical field of medical robots, in particular to a lung bronchus voxel model optimization method for monocular SLAM dynamic modeling. BACKGROUND

[0002] A lung intervention surgery robot diagnosis and treatment system constructs a cavity environment model and locates the position of the robot in the body at the same time, and then plans a path for the robot to perform surgery. A lung bronchus voxel map sequence is established based on monocular SLAM, and the surgery robot can plan a movement according to whether a voxel block in the map sequence can pass through. However, the estimated lung bronchus sparse point cloud map is inaccurate and contains many abnormal points, so that a large area on the map cannot be interpolated during the densification process, and finally a voxel map with surface pores is obtained. Due to the error of the constructed lung bronchus voxel map, the lung cavity model cannot accurately assist the lung intervention surgery robot when planning and navigating.

[0003] In theory, only one voxel block is needed to describe the inner wall of the lung bronchus cavity, however, the depth information of the cavity feature points estimated by monocular SLAM has errors, and the large area is not interpolated during the densification process of the sparse point cloud map, resulting in the phenomenon of redundancy and pores in the constructed lung cavity voxel map. The lung bronchus has a complex shape, and direct linear fusion interpolation of the entire map cannot accurately optimize the model. In order to solve this problem, a lung bronchus voxel model optimization method for monocular SLAM dynamic modeling is studied, and a radial projection interpolation strategy is proposed to optimize the lung bronchus voxel model. SUMMARY

[0004] The application aims to provide a lung bronchus voxel model optimization method for monocular SLAM dynamic modeling, which solves the problems of voxel redundancy and surface pores in the voxel model constructed by the lung intervention surgery robot in the lung dynamic bronchus environment.

[0005] The technical scheme of the application is as follows: a lung bronchus voxel model optimization method for monocular SLAM dynamic modeling, comprising the following steps:

[0006] Step 1, when the lung intervention surgery robot is inserted into the lung bronchus, it moves along the cavity center, constructs an environment model based on monocular SLAM, simultaneously locates the trajectory of the front end probe of the robot, and uses the robot trajectory as an approximate cavity center line;

[0007] Step 2, establish the tangent plane of the robot trajectory at a voxel unit interval, project the voxel block to the tangent plane to form an image describing the radial profile of the cavity, and the pixels in the image correspond to the voxel block in the lung bronchus voxel map through the projection model;

[0008] Step 3, using a radial line segment search algorithm to eliminate the radial redundant voxels of the lung bronchus model, performing graph interpolation based on curvature consistency to complete the filling of the model surface pores;

[0009] Step 4, according to the graph interpolation result based on curvature consistency, based on the mapping relationship between two-dimensional pixels and three-dimensional voxel blocks, the optimized two-dimensional grid is converted into a three-dimensional bronchus voxel map.

[0010] In the above technical solution, the pose of the lung intervention surgery robot in step 1 is represented in the form of three-dimensional coordinates and quaternions. The three-dimensional coordinates represent the translation component in the pose, and the quaternion represents the rotation component. The quaternion includes a real part and three imaginary parts, which are defined as follows:

[0011] q=q0+q1i+q2j+q3k (1)

[0012] In the above formula, q0 is the real part, q1, q2, and q3 are the imaginary parts, and i, j, and k are the components in the direction of the three rotation axes. The first frame in the robot trajectory is the reference frame in SLAM, which is the origin of the world coordinate system by default. The key frame poses after that are the translation and rotation relative to the first frame. The key frame pose represents the movement and orientation of the camera lens, using the translation vector t and the rotation matrix R world→cam respectively. world→cam The rotation matrix R

[0013]

[0014] The camera coordinates of a point in three-dimensional space are Pc, which are transformed into world coordinates P i by the translation vector t i cam→world and the rotation matrix R w , which can be represented by the following formula:

[0015] P w = R i cam→world P c +t i (3)

[0016] In the above formula, the rotation matrix i represents the i-th key frame in the trajectory. To simplify the calculation, the non-homogeneous coordinates of the space point are converted into homogeneous coordinates, and formula (3) can be rewritten as follows:

[0017]

[0018] In the above formula, and are the homogeneous coordinate forms of P w and P c , respectively, T i cam→world represents the i-th key frame transformation matrix, which is composed of a translation vector t i and a rotation matrix R i cam→world .

[0019] Thus, the lung intervention surgery robot trajectory is taken as the approximate cavity center line, and the transformation matrix T i from the camera coordinate system to the world coordinate system of a frame in the robot key frame trajectory at the corresponding time is obtained by using the world coordinate t i world→cam and the rotation matrix R i cam→world of the frame, so as to complete the coordinate transformation and lay a foundation for the subsequent projection model establishment and lung bronchus voxel model optimization.

[0020] In the above technical solution, the coordinate system of the tangent plane established in step 2 takes the top left corner vertex of the image as the origin O, the width as the X axis, and the height as the Y axis. Set ncol as the image width, nrow as the image height, and the pixel size as rs*rs, wherein rs is less than or equal to the edge length of the voxel in the cavity voxel map, so as to prevent multiple voxels from being projected into the same pixel. The Z axis of the camera coordinate system passes through the image center of the tangent plane, and the X axis and the Y axis of the camera coordinate system are projected on the image plane as x c and y c . The coordinates of a pixel in the image are (u, v), the camera coordinates of the pixel in the three-dimensional space are (x c , y c , z c ), and the camera coordinates of the corresponding voxel of the pixel can be obtained, and the calculation formula is as follows:

[0021]

[0022] In the above formula, Z ij represents the distance between the j-th tangent plane extracted in the camera coordinate system corresponding to the i-th key frame and the camera origin.

[0023] In the above technical solution, the image projected onto the voxel in the cross-section in step 3 can be regarded as a two-dimensional grid map, with pixels as occupied or vacant grids, corresponding to the voxel state. A radial line segment search method is used to traverse all grids in the cross-section image and query the occupied state of all grids. Boundary points (mincol, minrow), (maxcol, minrow), (maxcol, maxrow), and (mincol, maxrow) are set according to the center coordinate range of the occupied grids, and connected sequentially to obtain the boundary line. The center point 0 of the image is taken as a fixed point of the search line segment, and the point (mincol, minrow) on the boundary line is taken as the initial starting point of the line segment. Connecting point 0 and point (mincol, minrow) yields all occupied grids intersecting with this line segment. Then, the starting point of the search line segment on the boundary line is changed, and the search proceeds by moving one grid clockwise along the boundary line each time, until all occupied grids have been traversed.

[0024] In the above technical solution, step 3 eliminates radially redundant voxels in the pulmonary bronchus model. A radial line segment search algorithm is used to obtain the set of occupied gratings W on the cross-section, where occupied gratings on a search line segment are stored in a subset K. If the number of gratings stored in K is greater than 1, it indicates radial voxel redundancy in the cavity model corresponding to the occupied gratings along that line segment direction; in this case, multiple gratings on the line segment are merged into one. The distance d from the center coordinates of each occupied grating on the search line segment to point O is calculated. j And based on the occupancy probability w of each grid cell j Estimate the distance d0 from the merged raster to the center point O, and take d j The difference Δd between d0 and d0 j When Δd j When the value is minimized, the corresponding raster is the fused result. The radial redundancy voxels in the cavity are eliminated using the above data fusion method, at which point the number of radial raster grids in the voxel projection onto the cross-section is ≤1.

[0025] In the above technical solution, step 3 uses curvature consistency graphic interpolation to fill the pores on the cavity surface. A radial line segment search method is used to traverse the occupied grid, resulting in a set U =

[0026] {V1,V2,V3,...,V n}, where subset V i The grid stores the occupied grid cells on the i-th search line segment. The graph curve equation y is fitted using the coordinates of the occupied grid cells. The curvature at the center of each occupied grid cell is calculated based on y, and the mean curvature K is taken, as shown in the following formula:

[0027]

[0028] In the above formula, N represents the number of grid cells occupied in the graphic. denotes the first derivative of the curve at the i-th occupied grid, denotes the second derivative of the curve at the i-th occupied grid. Thus, the graph radius r is calculated as follows:

[0029]

[0030] When the number of occupied grids in the subset V i is 0, it indicates that there is a hole on the corresponding cavity wall at this position, which needs to be interpolated. In the subset V i , the grid at r distance from the center point O is added on the corresponding search line segment, and the point corresponding to the added grid is the interpolation point. Then, all the extracted sections are interpolated based on the curvature consistency interpolation method.

[0031] In the above technical solution, in step 4, the two-dimensional grid converted to the three-dimensional bronchial voxel map according to the mapping relationship between the two-dimensional plane pixel and the three-dimensional voxel unit in step 2 is optimized. The pixel coordinates in the section image are (u, v), and the camera coordinates in the three-dimensional space to which the pixel coordinates are mapped can be calculated according to formula (6) as (x c ,y c ,z c ). Then, according to the transformation relationship between the camera coordinate system and the world coordinate system in step 1, the point (x c ,y c ,z c ) is transformed to the world coordinate system according to formula (3), and finally the optimized lung bronchial voxel model is obtained.

[0032] The advantages of the present application are:

[0033] The present application provides a lung interventional surgery robot based on monocular SLAM dynamic modeling of lung bronchial voxel model optimization method, which estimates the key frame trajectory of the lung interventional surgery robot based on monocular SLAM instead of the cavity center line, then extracts the bronchial voxel model radial section perpendicular to the robot trajectory at intervals of the minimum voxel edge length, and establishes the projection model of the model in the radial direction of the bronchus. Then, based on the radial line segment search method, the occupied grids on the section are traversed, a curvature consistency graph interpolation method is proposed, the bronchial model radial redundant voxel elimination and surface hole filling are realized. Finally, based on the mapping relationship between the two-dimensional pixel and the three-dimensional voxel block, the section data is converted to the voxel map, and the optimized lung bronchial voxel model is obtained. The present application can solve the problems of voxel redundancy and surface hole in the voxel model. BRIEF DESCRIPTION OF DRAWINGS

[0034] The present application will be further described below in conjunction with the drawings and examples:

[0035] Figure 1 Flow chart of the method of embodiment one of the present application.

[0036] Figure 2 Extraction diagram of the radial section of the bronchial model according to the present application.

[0037] Figure 3 Diagram of the section image coordinate system constructed by the method according to the present application.

[0038] Figure 4 Diagram of the section image of the voxel model constructed by the method according to the present application.

[0039] Figure 5 Diagram of the radial line segment search method according to the present application.

[0040] Figure 6 Radial section graph interpolation result diagram according to the present application. DETAILED DESCRIPTION

[0041] In order to achieve the purpose of the present application, as shown in the figure, in one embodiment of the present application, a lung bronchial voxel model optimization method for monocular SLAM dynamic modeling includes the following steps: Figure 1

[0042] In step 1, the pose of the lung interventional surgery robot is represented in the form of three-dimensional coordinates and quaternions. The three-dimensional coordinates represent the translation component in the pose, and the quaternion represents the rotation component. The quaternion includes a real part and three imaginary parts, which are defined as follows:

[0043] q = q0 + q1i + q2j + q3k (1)

[0044] In the above formula, q0 is the real part, q1, q2, and q3 are the imaginary parts, and i, j, and k are the components in the direction of the three rotation axes. The first frame in the robot trajectory is the reference frame in SLAM, which is the origin of the world coordinate system by default. The key frame poses after that are the translation and rotation relative to the first frame. The key frame pose represents the movement and orientation of the camera lens, using the translation vector t and the rotation matrix R world→cam Description, R world→cam which is calculated by the following formula:

[0045]

[0046] The camera coordinates of a point in three-dimensional space are P c , which is transformed into the world coordinates P i by the translation vector t i cam→world and the rotation matrix R w , which can be expressed by the following formula:

[0047] P w = R​i cam→world P c +t i (3)

[0048] In the above formula, the rotation matrix i represents the i-th key frame in the trajectory. To simplify the calculation, the non-homogeneous coordinates of the spatial point are converted into homogeneous coordinates, and formula (3) can be rewritten as follows:

[0049]

[0050]

[0051] In the above formula, and are the homogeneous coordinate forms of P w and P c , respectively, T i cam→world represents the i-th key frame transformation matrix, which is composed of the translation vector t i and the rotation matrix R i cam→world .

[0052] Thus, the lung intervention surgery robot trajectory is taken as the approximate lumen centerline, and the transformation matrix T i cam→world from the camera coordinate system to the world coordinate system of a frame in the robot key frame trajectory at the corresponding time is obtained by the world coordinate t i and the rotation matrix R i world→cam of the frame, the coordinate transformation is completed, and the basis for the subsequent projection model establishment and lung bronchus voxel model optimization is established.

[0053] In step 2, the tangent plane of the robot trajectory is established at a voxel unit interval, as shown in Figure 2 . The voxel block is projected onto the tangent plane to form an image describing the radial profile of the lumen, and the pixels in the image correspond to the voxel blocks in the lung bronchus voxel map through the projection model. The coordinate system of the tangent plane is shown in Figure 3 , with the top left corner of the image as the origin O, the width as the X axis, and the height as the Y axis. Set ncol as the image width, nrow as the image height, and the pixel size as rs*rs, where rs is less than or equal to the edge length of the voxel in the lumen voxel map to prevent multiple voxels from being projected into the same pixel. The camera coordinate system Z axis passes through the image center of the tangent plane, and the camera coordinate system X axis and Y axis are projected onto the image plane as x c and y c . The coordinates of a pixel in the image are (u, v), and its mapping to the camera coordinates in three-dimensional space can be calculated by formula (6) as (x c , y c , z c ). The coordinates of the corresponding voxel of this pixel in the camera coordinate system can be obtained, and the calculation formula is as follows:

[0054]

[0055] In the above formula, Z ij This represents the distance between the j-th slice extracted in the camera coordinate system corresponding to the i-th keyframe and the camera origin.

[0056] In step 3 Figure 4 This is a projection image of a bronchial voxel map of the lungs onto a radial section, approximating an "O" shape. The boxes in the image indicate voxel redundancy in the cavity walls, while the circles show larger pores on the cavity map.

[0057] The image projected onto a voxel in a cross-section can be viewed as a two-dimensional raster map, with pixels representing occupied or vacant raster cells, corresponding to the voxel's state. A radial line segment search method is used to traverse all raster cells in the cross-section image and query the occupancy state of each cell. For example... Figure 5 As shown, the black circles represent occupied grids in the image plane. Boundary points (mincol, minrow), (maxcol, minrow), (maxcol, maxrow), and (mincol, maxrow) are set according to the center coordinates of the occupied grids, and these are connected sequentially to obtain the boundary line, as shown in the gray area of ​​the figure. Using the center point O of the image as a fixed point for the search line segment, and the point (mincol, minrow) on the boundary line as the initial starting point of the line segment, connecting point O and point (mincol, minrow) retrieves all occupied grids intersecting with this line segment. Then, the starting point of the search line segment on the boundary line is changed, and the search proceeds by moving one grid clockwise along the direction of the dotted line in the figure, i.e., the boundary line, until all occupied grids have been traversed.

[0058] To eliminate radial redundancy voxels in the pulmonary bronchial model, a radial segment search algorithm is used to obtain the set of occupied gratings W on the cross-section. Occupied gratings on a search segment are stored in a subset K. If the number of gratings stored in K is greater than 1, it indicates radial voxel redundancy in the cavity model corresponding to the occupied gratings along that segment direction; in this case, multiple gratings on the segment are merged into one. The distance d from the center coordinates of each occupied grating on the search segment to point O is calculated. j And based on the occupancy probability w of each grid cell j Estimate the distance d0 from the merged raster to the center point O, and take d j The difference Δd between d0 and d0 j When Δd j When the value is minimized, the corresponding raster is the fused result. The radial redundancy voxels in the cavity are eliminated using the above data fusion method, at which point the number of radial raster grids in the voxel projection onto the cross-section is ≤1.

[0059] Based on the curvature consistency graphic interpolation, the cavity surface pore filling is realized. The radial line segment search method is used to traverse the occupancy grid, and the set U = {V1, V2, V3,..., V n} is obtained, wherein the sub-set V i stores the occupancy grid on the i-th search line segment. The occupancy grid coordinates are used to fit the graphic curve equation y, the curvature at the center of each occupancy grid is calculated according to y, and the average curvature K is taken, as follows:

[0060]

[0061] In the above formula, N represents the number of occupancy grids in the graph, represents the first derivative of the curve at the i-th occupancy grid, represents the second derivative of the curve at the i-th occupancy grid. Thus, the radius r of the graph is calculated as follows:

[0062]

[0063] When the number of occupancy grids in the sub-set Vi is 0, it indicates that there is a pore on the cavity wall corresponding to this place, and interpolation is needed. In the sub-set V i with 0 occupancy grids, a grid is added at r from the center point O on the corresponding search line segment, and the point corresponding to the added grid is the interpolation point. Figure 6 The graphic interpolation result is shown in the figure, and the hollow part of the graph is filled by interpolation to complete the cavity pore filling. Then, all the extracted sections are interpolated based on the curvature consistency interpolation method.

[0064] In step 4, according to the mapping relationship between the two-dimensional plane pixels and the three-dimensional voxel units in step 2, the optimized two-dimensional grid is converted into a three-dimensional bronchial voxel map. The pixel coordinates in the section image are (u, v), and the camera coordinates in the three-dimensional space can be calculated by formula (6) as (x c , y c , z c ). Then, according to the transformation relationship between the camera coordinate system and the world coordinate system in step 1, the point (x c , y c , z c ) is transformed to the world coordinate system by formula (3), and finally the optimized lung bronchial voxel model is obtained.

[0065] Of course, the above embodiments only illustrate the technical concept and characteristics of the present application, and the purpose is to enable those skilled in the art to understand the content of the present application and implement it, and cannot limit the protection scope of the present application. Any modification made according to the spirit and essence of the main technical solution of the present application should be covered within the protection scope of the present application.

Claims

1. A method for optimizing lung bronchial voxel models for dynamic modeling in monocular SLAM, characterized in that, Includes the following steps: Step 1: When the pulmonary interventional surgery robot is inserted into the bronchus of the lung, it moves along the center of the cavity. While constructing an environmental model based on monocular SLAM, the trajectory of the robot's front probe is simultaneously located, and the robot trajectory is used as an approximate center line of the cavity. Step 2: Using voxel units as intervals, establish a tangent plane for the robot trajectory, project the voxel blocks onto the tangent plane to form an image describing the radial contour of the cavity. The pixels in the image correspond to the voxel blocks in the lung bronchus voxel map through the projection model. Step 3: Use the radial segment search algorithm to eliminate radial redundant voxels in the lung bronchus model, and perform graphic interpolation based on curvature consistency to fill the pores on the model surface. Step 4: Based on the curvature consistency graphic interpolation results, and based on the mapping relationship between two-dimensional planar pixels and three-dimensional voxel blocks, convert the optimized two-dimensional raster into a three-dimensional bronchial voxel map. In step 3, the image of the voxel projection in the cross section is used as a two-dimensional grid map, and the pixels are used as occupied or empty grids, corresponding to the voxel state; based on the radial line segment search method, all grids of the cross section image are traversed to query the occupied state of all grids. The method of traversing all grid cells of the cross-sectional image based on radial line segment search and querying the occupancy status of all grid cells includes: Set boundary points (mincol, minrow), (maxcol, minrow), (maxcol, maxrow), and (mincol, maxrow) based on the center coordinate range of the occupied grid, and connect them sequentially to obtain the boundary line. Take the center point O of the graphic as a fixed point of the search line segment, and the point (mincol, minrow) on the boundary line as the initial starting point of the line segment. Connect point O and point (mincol, minrow) to obtain all occupied grids that intersect with the line segment. Change the starting point of the search line segment on the boundary line, and move one grid clockwise along the boundary line each time to search, until all occupied grids have been traversed. In step 3, eliminating radial redundancy voxels in the lung bronchial model includes: The radial segment search algorithm is used to obtain the set W of occupied grids on the cross-section, where the occupied grids on a search segment are stored in a subset K. If the number of grids stored in K is greater than 1, it indicates that the cavity model corresponding to the occupied grid in the direction of the segment is radially redundant, and the multiple grids on the segment are merged into one. The distance d from the center coordinates of each occupied grid on the search segment to point O is calculated. j And based on the occupancy probability w of each grid cell j Estimate the distance d0 from the merged raster to the center point O, and take d j The difference Δd between d0 and d0 j When Δd j When it is at its minimum, the corresponding grid is the result of fusion; In step 3, the step of performing graphic interpolation based on curvature consistency to fill the pores on the cavity surface includes: Using the radial line segment search method to traverse the occupied grid, we obtain the set U = {V1, V2, V3, ..., V}. n }, where subset V i The grid stores the occupied grid cells on the i-th search line segment; the graph curve equation y is fitted using the coordinates of the occupied grid cells, and the curvature at the center of each occupied grid cell is calculated based on y, and the mean curvature K is taken, as shown in the following formula: In the above formula, N represents the number of grid cells occupied in the graphic. This represents the first derivative of the curve at the i-th occupied grid cell. Let represent the second derivative of the curve at the i-th occupied grid cell; thus, the radius r of the graphic is calculated as follows: In the subset V that occupies 0 grid cells i On the corresponding search line segment, a grid is added at a distance r from the center point O. The point corresponding to this added grid is the interpolation point. All extracted cross-sections are interpolated based on the curvature consistency interpolation method.

2. The method for optimizing lung bronchial voxel models for monocular SLAM dynamic modeling according to claim 1, characterized in that: In step 1, the trajectory of the pulmonary interventional surgery robot is used as an approximate cavity centerline, and the world coordinates t of a certain frame in the robot's keyframe trajectory at the corresponding time are used. i and rotation matrix R iworld→cam The transformation matrix T from the camera coordinate system to the world coordinate system of this frame is obtained. icam→world Complete the coordinate transformation.

3. The method for optimizing lung bronchial voxel models for monocular SLAM dynamic modeling according to claim 2, characterized in that: In step 2, a coordinate system for the tangent plane is established, with the top-left corner of the image as the origin O, the width as the X-axis, and the height as the Y-axis. ncol is set as the image width, nrow as the image height, and the pixel size as rs*rs, where rs is less than or equal to the side length of a voxel in the cavity voxel map. The camera coordinate system's Z-axis passes through the image center of the tangent plane, and the projections of the camera coordinate system's X and Y axes onto the image plane are x... c and y c The coordinates of a pixel in the image are (u, v), and the camera coordinates mapped to this pixel in 3D space are (x, v). c ,y c ,z c The value of a pixel can be obtained by using the coordinates of its corresponding voxel in the camera coordinate system. The calculation formula is as follows: In the above formula, Z ij This represents the distance between the j-th slice extracted in the camera coordinate system corresponding to the i-th keyframe and the camera origin.

4. The method for optimizing lung bronchial voxel models for monocular SLAM dynamic modeling according to claim 3, characterized in that: Step 4 specifically includes: Let the pixel coordinates in the cross-sectional image be (u, v), and calculate the camera coordinates in three-dimensional space mapped from the pixel coordinates to (x, v) using equation (1). c ,y c ,z c Based on the transformation relationship between the camera coordinate system and the world coordinate system in step 1, the point (x... c ,y c ,z c The model is then transformed to the world coordinate system, resulting in the optimized lung bronchus voxel model.

Images