An intelligent analysis and surgery planning auxiliary system for thoracic surgery images

Abstract

The present application provides a kind of intelligent analysis and surgery planning auxiliary system of thoracic surgery image, obtains patient chest tomography image data, separates trachea-bronchial tree, blood vessel and lesion, constructs three-dimensional anatomical model;Blood vessel and bronchial wall local thickness, CT value determine anatomical barrier stiffness index, combined with blood vessel diameter and lung functional density distribution, give risk generation value to spatial voxel, construct operation risk potential field containing stiffness index positive correlation penalty term;Through the irregularity of lesion surface, determine the multi-point target area, search the initial path of body surface puncture point to target area in risk potential field, and the path iteration step and direction are controlled by the local gradient tensor of risk potential field;Calculate the local curvature of initial path, generate variable boundary safety channel according to the risk potential field gradient and curvature around path, and the radial width of channel is inversely proportional to risk potential field gradient and proportional to local curvature, and the initial path is optimized smoothly with channel boundary as constraint, to obtain surgery planning path.

Description

Technical Field

[0001] This application belongs to the field of path planning, and in particular relates to an intelligent analysis and surgical planning assistance system for thoracic surgical images. Background Technology

[0002] Manual planning based on 2D tomographic images or 3D reconstruction models is highly subjective, time-consuming, and makes it difficult to intuitively assess puncture risks within the anatomical structure. Computer-aided surgical planning technology reconstructs a 3D model by segmenting patient images and then plans the geometrically shortest path from the body surface to the lesion. Current technologies simplify anatomical structures into obstacles to be avoided in risk assessment models, failing to utilize the differences in characteristics of different tissues and structures. Path search strategies are performed on a point-by-point basis, with the target point being the geometric center of the lesion, failing to fully consider the size and shape of the surgical instruments themselves, and ignoring the accessibility and therapeutic value of different areas on the lesion surface. Regarding the presentation and optimization of planning results, most methods only provide an idealized path, failing to provide a safety boundary reference for intraoperative operations, and high-curvature inflection points in the path may not conform to the motion constraints of the instruments. Simple smoothing of the path, lacking safety boundary constraints, may cause the optimized path to deviate into insufficiently assessed risk areas. Therefore, there is an urgent need for a planning method that can construct a risk field and generate a safe passage that balances safety and operability. Summary of the Invention

[0003] This invention proposes an intelligent image analysis and surgical planning assistance system for thoracic surgery to address the problem that existing technologies cannot construct a risk field and are difficult to generate planning methods for safe passages that balance safety and operability. The system includes: The first construction module is used to acquire the patient's chest tomographic image data, segment the tracheobronchial tree, blood vessels and lesions, and construct a three-dimensional anatomical model. The second construction module is used to calculate the local thickness and CT value of the blood vessel and bronchial wall based on the three-dimensional anatomical model to determine the stiffness index of the anatomical barrier, and combine the blood vessel diameter and lung lobe functional density distribution to construct a surgical risk potential field by assigning risk value to the spatial voxels. The search module is used to calculate the surface irregularity of the lesion to determine the multi-point target area on the lesion surface; using the geometric model of the surgical instrument as the search unit, it searches for the initial path from the puncture point on the body surface to any point in the multi-point target area in the risk potential field; the iteration step size and direction of the search unit are controlled by the local gradient tensor of the risk potential field. An optimization module is used to calculate the local curvature of the initial path; and based on the risk potential gradient around the initial path and the local curvature, generate a variable boundary safe passage around the initial path; and use the boundary of the safe passage as a constraint to perform smooth optimization on the initial path to obtain the planned path.

[0004] This invention constructs a surgical risk potential field by representing the stiffness index of blood vessel and bronchial walls and integrating information on blood vessel diameter and lung lobe functional density. In path search, the geometric model of surgical instruments is used as the search unit, searching multiple target areas on the lesion surface. The search process is guided by the local gradient of the risk field, ensuring that the planned initial path closely matches the actual operation. A safe passage is generated around the initial path, with its width adjusted according to the surrounding risk gradient and path curvature, providing the surgeon with a safe operating range. Using this safe passage as a rigid constraint, the path is optimized to ensure that it does not deviate from high-risk areas after smoothing, improving the safety and clinical operability of surgical planning. Attached Figure Description

[0005] Figure 1 This is a schematic diagram of a three-dimensional anatomical model; Figure 2 A schematic diagram of the surgical risk potential field; Figure 3 This is a schematic diagram of path optimization and safe passage. Detailed Implementation

[0006] The features and exemplary embodiments of various aspects of this application will be described in detail below. To make the objectives, technical solutions, and advantages of this application clearer, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are only intended to explain this application and not to limit it. For those skilled in the art, this application can be implemented without some of these specific details. The following description of the embodiments is merely to provide a better understanding of this application by illustrating examples.

[0007] It should be noted that, in this document, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising..." does not exclude the presence of additional identical elements in the process, method, article, or apparatus that includes said element.

[0008] In the first embodiment, the present invention proposes an intelligent analysis and surgical planning assistance system for thoracic surgical images, comprising: The first construction module is used to acquire the patient's chest tomographic image data, segment the tracheobronchial tree, blood vessels and lesions, and construct a three-dimensional anatomical model. DICOM format chest CT images of the patient were acquired using medical equipment. A region-growing algorithm was employed to segment the complete trachea-bronchial tree structure by setting seed points within the trachea. A Hessian matrix-based vascular enhancement filtering algorithm was used to enhance the vascular structure, and a thresholding method was then used to extract the complete pulmonary vascular network. For lesions, a pre-trained U-Net deep learning network was used for segmentation. The binary masks of each segmented tissue and organ were then processed using the Marching Cubes algorithm to generate their respective triangular mesh surface models, which together constituted a three-dimensional anatomical model, such as... Figure 1 .

[0009] The second construction module is used to calculate the local thickness and CT value of the blood vessel and bronchial wall based on the three-dimensional anatomical model to determine the stiffness index of the anatomical barrier, and combine the blood vessel diameter and lung lobe functional density distribution to construct a surgical risk potential field by assigning risk value to the spatial voxels. The local thickness of blood vessel and bronchial walls is obtained by calculating the Euclidean distance from each point on the surface of the 3D model to the corresponding point on the inner or outer wall of the model. Simultaneously, the Heinz unit values ​​at the corresponding barrier locations are sampled from the original CT images as CT values. The stiffness index is defined as the weighted sum of the local thickness and the CT value. The blood vessel diameter is calculated by fitting the maximum inscribed sphere at each point on the centerline of the blood vessel model. Lung lobe functional density is determined by analyzing the CT value distribution of the lung parenchyma region; regions with lower CT values ​​represent higher air content, better function, and correspondingly higher risk assignments. The risk cost of any voxel in space is set as a comprehensive function, the value of which is inversely proportional to the distance from the voxel to the blood vessel and bronchus, directly proportional to the diameter of adjacent blood vessels, and included in the lung lobe functional density value. When the path planning algorithm determines that crossing the blood vessel or bronchial wall is necessary, a penalty value proportional to the barrier stiffness index is added to the total cost.

[0010] In an optional embodiment, the step of calculating the local thickness and CT values ​​of the vascular and bronchial walls based on the three-dimensional anatomical model to determine the stiffness index of the anatomical barrier includes: For each surface point p on the wall of a blood vessel or bronchus in the three-dimensional anatomical model, the local thickness T(p) is measured using the normal ray projection method, and the CT value H(p) corresponding to the point is extracted. The thickness and CT values ​​are normalized to obtain the normalized thickness. and normalized CT value ; The stiffness exponent S(p) of the point is calculated using the following formula: in, and These are the preset weighting coefficients.

[0011] Specifically, a three-dimensional mesh model of the blood vessels and bronchial walls is constructed based on the patient's CT image data. For any vertex p on the model surface, two virtual rays are emitted inward and outward along the normal direction of that point until the rays penetrate the wall on the other side. The straight-line distance between the two intersection points is defined as the local thickness T(p). For example, if the thickness of a point p on a bronchial wall is measured to be 2.5 mm, and the original CT data is consulted, the CT value at point p is found to be 120 Henle units.

[0012] Assuming the lung wall thickness ranges from 1 mm to 8 mm and the CT value ranges from -50 Henlein units to 300 Henlein units in the entire lung anatomical model, then the normalized thickness at point p is... Approximately equal to 0.214, normalized CT value Approximately 0.486. Weighting coefficients are set to balance the contributions of thickness and density to stiffness, for example, by setting... It is 0.6. If the value is 0.4, then the stiffness exponent S(p) at point p is calculated to be approximately 0.323. This process is performed on all points on the barrier surface to obtain a complete stiffness distribution map.

[0013] In an optional embodiment, the functional density of the lung lobe is quantitatively calculated based on the distribution of CT values ​​in the lung parenchyma region of the patient's chest CT image: First, the complete lung parenchyma region is extracted using image segmentation technology, and non-lung parenchyma structures such as blood vessels, bronchi, and lesions are removed; the range of CT values ​​for all voxels in this region is statistically analyzed, and the min-max normalization method is used to map the CT values ​​of each voxel to dimensionless values ​​in the interval [0,1]; since a lower CT value represents a higher air content and better gas exchange function in the lung tissue, the normalized value is directly used as the functional density of the corresponding voxel, F(v). The better the function of the region, the closer F(v) is to 1. In the subsequent risk potential field construction, this value will be included as a positive weight term in the voxel risk cost value to achieve key protection of functionally important lung tissues.

[0014] To represent the surgical space as a three-dimensional risk map, in an optional embodiment, the process of constructing a surgical risk potential field by assigning risk values ​​to spatial voxels includes:

[0015] For any voxel v in space, the risk cost R(v) of the voxel is determined by the following formula: in, The stiffness index is the stiffness index of the anatomical barrier closest to voxel v. The diameter of the blood vessel closest to voxel v. This represents the functional density of the lung lobe at the location of voxel v. , These are the dimensionless values ​​of each term after min-max normalization. , , These are positive weighting coefficients.

[0016] For each three-dimensional voxel v within the lung CT scan coverage area, the voxel risk value R(v) is a weighted sum of three parts. The first part is the stiffness index. Calculate the shortest distance from voxel v to the surfaces of all blood vessels and bronchi, and adopt the stiffness index S(p) of the nearest surface point as the voxel's... For example, if the nearest structure to a voxel v is a section of bronchial wall, and the stiffness index of the nearest point on that wall is 0.4, then the voxel's... That is, 0.4. The second part is the vessel diameter, which calculates the shortest distance from voxel v to the centerline of all vessels, and adopts the diameter of the nearest vessel as the mean. For example, if voxel v is closest to an artery with a diameter of 4 mm, then... The value is 4. The third part is the functional density of the lung lobe, F(v). This data is usually obtained from functional imaging or preset according to anatomical location. For example, if voxel v is located in the functionally important left lower lobe, the functional density F(v) can be set to 0.9, while the voxel located in the less functional area may be 0.3. The values ​​are obtained after min-max normalization. , By setting weights, such as Calculate the total risk value R(v) for voxel v. Perform this calculation on all voxels to generate the complete surgical risk potential field, such as... Figure 2 .

[0017] In one optional embodiment, the functional density of the lung lobe is quantitatively calculated based on the distribution of CT values ​​in the lung parenchyma region of the patient's chest CT images: First, the complete lung parenchyma region is extracted using image segmentation technology, and non-parenchymal structures such as blood vessels, bronchi, and lesions are removed; then, the CT value range of all voxels within this region is statistically analyzed, and the min-max normalization method is used to map the CT value of each voxel to a dimensionless value in the interval [0,1].

[0018] The search module is used to calculate the surface irregularity of the lesion to determine the multi-point target area on the lesion surface; using the geometric model of the surgical instrument as the search unit, it searches for the initial path from the puncture point on the body surface to any point in the multi-point target area in the risk potential field; the iteration step size and direction of the search unit are controlled by the local gradient tensor of the risk potential field.

[0019] The surface irregularity of the lesion is represented by calculating the average curvature of each vertex on the 3D surface model. The set of vertices with curvature greater than a preset threshold is defined as the multi-point target region. Surgical instruments, such as puncture needles, are abstracted as cylindrical models with specific lengths and radii, with the central axis of the cylinder as the search target. An improved A* search algorithm is used, starting the search from the puncture point specified by the doctor on the body surface model. In each iteration, the algorithm calculates the risk potential gradient tensor at the current search unit's location and prioritizes exploring along the eigenvector direction corresponding to the smallest eigenvalue of the gradient tensor. Simultaneously, the iteration step size is adjusted according to the gradient magnitude; the larger the magnitude, the smaller the step size. The search stops when any part of the search unit first contacts any vertex in the multi-point target region, and the search trajectory is recorded as the initial path.

[0020] To identify the most likely malignant areas on the lesion surface for biopsy or ablation, in an optional embodiment, calculating the surface irregularity of the lesion to determine multi-point target areas on the lesion surface includes: Calculate the average curvature of each vertex on the surface mesh of the 3D model of the lesion; The set of vertices with an average curvature greater than a preset threshold is defined as candidate target points; The candidate target points are clustered into 5 clusters using the K-means clustering algorithm, and the centroid of each cluster is used as a target point in the multi-point target region.

[0021] Specifically, a three-dimensional surface mesh model is constructed from the lesions segmented from the CT image. This model consists of thousands of small triangular patches and vertices. For each vertex, the average curvature of that point is calculated based on its own position and the positions of its adjacent vertices. This value represents the degree of surface roughness at that point. For example, a region with a smooth, spherical surface will have a lower average curvature value at its vertices, such as 0.05; while a region with sharp protrusions or burrs will have a higher average curvature value at its vertices, such as reaching 0.8.

[0022] A curvature threshold, such as 0.6, is set to filter out all vertices with an average curvature greater than 0.6, forming a set of candidate target points. These points represent the most irregular locations on the lesion surface. Assuming 300 candidate points are selected, to avoid overly dense clustering, the K-means clustering algorithm is used to process the 3D coordinates of these 300 points. The number of clusters, K, is set to 5, dividing the 300 points into 5 spatially clustered groups. The average 3D coordinates of all points within each cluster are calculated to obtain the centroid of that cluster. These 5 centroids are then identified as the multi-point target region in the surgical planning.

[0023] To achieve intelligent obstacle avoidance during path search, in an optional embodiment, the iteration step size and direction of the search unit are controlled by the local gradient tensor of the risk potential field, including: At the current position of the search unit, calculate the gradient tensor G of the risk potential field at that point; Perform eigenvalue decomposition on the gradient tensor G to obtain eigenvalues. and the corresponding feature vector ; will be with the minimum eigenvalue Corresponding feature vector The priority iteration direction for the search unit; the iteration step size L is set by the following formula: in, The base step size is k, and the adjustment coefficient is k.

[0024] When the search unit of a pathfinding algorithm is located at a point in space, it needs to determine the direction and distance of its next move. To do this, the gradient tensor G of the risk potential field at that point is calculated. This is a 3×3 matrix representing the rate of change and curvature of the risk value in each direction. For example, in a flat, low-risk region, the elements of G are close to zero; while next to a large blood vessel, the risk value increases along the direction perpendicular to the vessel, and the corresponding element values ​​of G will be very large.

[0025] Eigenvalue decomposition is performed on the gradient tensor G, yielding three eigenvalues ​​and their corresponding eigenvectors. These represent the three mutually orthogonal directions of risk change: the fastest, the second fastest, and the slowest, corresponding to the eigenvalues. This indicates the degree of drastic change in risk in the stated direction. Because... Minimum, corresponding direction This is the direction with the most gradual increase in risk, therefore it is chosen as the safest next move. Meanwhile, the iteration step size L is determined by the largest eigenvalue. Regulation. Assume a base step size. The value is 2mm, and the adjustment coefficient k is 0.5. If the current position is... A value of 3.0 indicates a strong risk source nearby, so the step size L is shortened to 0.8 mm. Conversely, if... It is very small, and the step size is close to the basic step size, enabling fast search.

[0026] An optimization module is used to calculate the local curvature of the initial path; and based on the risk potential gradient around the initial path and the local curvature, generate a variable boundary safe passage around the initial path; and use the boundary of the safe passage as a constraint to perform smooth optimization on the initial path to obtain the planned path.

[0027] The initial path is discretized into a series of path points. For each path point, a close circle is defined by the point and its two adjacent path points. The reciprocal of the radius of this close circle is the local curvature of that point. A safety radius is defined at each path point, which is proportional to the local curvature and inversely proportional to the gradient magnitude of the surrounding risk potential field. The safety radii of all path points are connected to form a variable boundary safety channel. Subsequently, the initial path is represented as a discrete set of path points P = {p1, p2, ..., pn}, and a target energy function containing smoothing and boundary constraint terms is constructed. The gradient descent method is used to iteratively adjust the positions of path points. A smoothing term pushes path points towards lower curvature to reduce tortuosity, while boundary constraints impose a penalty on points exceeding the channel to ensure safety. This process continues until the energy function converges. Finally, the optimized path points are fitted to generate a smooth B-spline curve, which serves as the planned path. Figure 3 As shown.

[0028] In an optional embodiment, generating a variable-boundary safe passage around the initial path based on the risk potential gradient around the initial path and the local curvature includes: For each point p on the initial path, calculate the local curvature at that point. and the mean gradient of the risk potential field in the local region surrounding the point ; The radial width W(p) of the channel at the point is calculated using the following formula: in, Based on the width, and It is a positive proportionality coefficient; A circular cross section is constructed on the normal plane with the radial width W(p) of each point on the path, and all the circular cross sections are connected to form the safety passage.

[0029] Specifically, iterate through every point p on the initial path. At each point p, calculate the curvature of the path at that point, i.e., the local curvature. For example, the curvature of a point p1 on a straight section of a path. The curvature is close to 0, but at point p2, which is at a sharp bend, the curvature is... It could be 0.9. Simultaneously, calculate the average value of the risk potential gradient within a small volume surrounding point p. This reflects the degree of risk change near that point. For example, p1 is located in open lung parenchyma, and the mean risk gradient may be only 0.1; while p2 is located close to an aorta, and the mean risk gradient may be as high as 5.0.

[0030] Assuming base width 4mm, adjustment coefficient It is 0.8. The width is 0.3. At point p1, the calculated width W(p1) is approximately 3.88 mm. At point p2, the calculated width W(p2) is approximately 2.75 mm. The passage will be appropriately widened at path bends to facilitate the passage of equipment, while the passage will narrow near hazardous structures to provide strict protection. A disk with a radius of this width is constructed on the normal plane at each point on the path, and all the disks are connected to form a safety pipe with varying width.

[0031] To optimize an initial path that may contain unnecessary twists and turns into a smooth path that is strictly within the safety passage, in an optional embodiment, the initial path is smoothed and optimized using the boundary of the safety passage as a constraint to obtain the planned path, including: The initial path is represented as a series of discrete path points. ; Construct the target energy function: Among them, the smoothing term The integral of the total curvature of the path, and the boundary constraint terms. The penalty function is the distance from the path point to the safety passage boundary, and the value of the function increases when the path point exceeds the passage boundary; The gradient descent method is used to iteratively adjust the positions of the path points until the target energy function E(P) converges to the minimum value, thus obtaining the planned path.

[0032] Specifically, the initial path consists of n discrete points. The core of the optimization is to minimize an energy function E(P), which contains two terms: a smoothing term and a... The smoothness of a path is measured, for example, by calculating and summing the squares of the curvature at each point on the path. A winding path will have high smoothness energy, while a path that is close to a straight line will have low energy.

[0033] Boundary constraint terms Ensure the path remains within the safe passage. For each point on the path. Calculate its distance to the nearest safe passage boundary. The energy value of this term is zero as long as the point is within the passage. Once a point... In the iterations, moving a point outside the channel boundary, even by just 0.1 mm, incurs a very large energy penalty. The optimization process employs gradient descent. In each iteration, the gradient of the energy function with respect to the coordinates of each path point is calculated, causing each point to move a small step in the direction that minimizes the total energy. For example, a point located at a sharp bend in the path will be pushed outward by the smoothing term, straightening the path, while a point too close to the channel boundary will be pushed back to the center of the channel by the repulsive force of the boundary term. This iterative process is repeated until the total energy no longer decreases; the resulting path point sequence P represents the planned smooth and safe path.

[0034] It should be clarified that this application is not limited to the specific configurations and processes described above and shown in the figures. For the sake of brevity, detailed descriptions of known methods are omitted here. In the above embodiments, several specific steps are described and shown as examples. However, the method process of this application is not limited to the specific steps described and shown. Those skilled in the art can make various changes, modifications, and additions, or change the order of steps, after understanding the spirit of this application.

[0035] The functional modules shown in the above-described block diagram can be implemented as hardware, software, firmware, or a combination thereof. When implemented in hardware, they can be, for example, electronic circuits, application-specific integrated circuits (ASICs), appropriate firmware, plug-ins, function cards, etc. When implemented in software, the elements of this application are programs or code segments used to perform the required tasks. Programs or code segments can be stored on a machine-readable medium or transmitted over a transmission medium or communication link via data signals carried on a carrier wave. "Machine-readable medium" can include any medium capable of storing or transmitting information. Examples of machine-readable media include electronic circuits, semiconductor memory devices, ROM, flash memory, erasable ROM (EROM), floppy disks, CD-ROMs, optical disks, hard disks, fiber optic media, radio frequency (RF) links, etc. Code segments can be downloaded via computer networks such as the Internet, intranets, etc.

[0036] It should also be noted that the exemplary embodiments mentioned in this application describe methods or systems based on a series of steps or apparatus. However, this application is not limited to the order of the above steps; that is, the steps can be performed in the order mentioned in the embodiments, or in a different order, or several steps can be performed simultaneously.

[0037] The aspects of this application have been described above with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It should be understood that each block in the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing apparatus to produce a machine such that these instructions, executable via the processor of the computer or other programmable data processing apparatus, enable the implementation of the functions / actions specified in one or more blocks of the flowchart illustrations and / or block diagrams. Such a processor can be, but is not limited to, a general-purpose processor, a special-purpose processor, a special application processor, or a field-programmable logic circuit. It is also understood that each block in the block diagrams and / or flowcharts, and combinations of blocks in the block diagrams and / or flowcharts, can also be implemented by dedicated hardware performing the specified functions or actions, or can be implemented by a combination of dedicated hardware and computer instructions.

[0038] The above description is merely a specific implementation of this application. Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems, modules, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here. It should be understood that the protection scope of this application is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in this application, and these modifications or substitutions should all be covered within the protection scope of this application.

Claims

1. A thoracic surgery image intelligent analysis and surgical planning assistance system, characterized in that, Includes the following modules: The first construction module is used to acquire the patient's chest tomographic image data, segment the tracheobronchial tree, blood vessels and lesions, and construct a three-dimensional anatomical model. The second construction module is used to calculate the local thickness and CT value of the blood vessel and bronchial wall based on the three-dimensional anatomical model to determine the stiffness index of the anatomical barrier, and combine the blood vessel diameter and lung lobe functional density distribution to construct a surgical risk potential field by assigning risk value to the spatial voxels. The search module is used to calculate the surface irregularity of the lesion and determine the multi-point target area on the lesion surface; Using the geometric model of the surgical instrument as the search unit, an initial path from the puncture point on the body surface to any point in the multi-point target region is searched in the risk potential field. The iteration step size and direction of the search unit are controlled by the local gradient tensor of the risk potential field. An optimization module is used to calculate the local curvature of the initial path; and based on the risk potential gradient around the initial path and the local curvature, to generate a variable boundary safe passage around the initial path. Using the boundaries of the safety passage as constraints, the initial path is smoothly optimized to obtain the planned path.

2. The system according to claim 1, characterized in that, The calculation of the local thickness and CT values ​​of the vascular and bronchial walls based on the three-dimensional anatomical model to determine the stiffness index of the anatomical barrier includes: For each surface point p on the wall of a blood vessel or bronchus in the three-dimensional anatomical model, the local thickness T(p) is measured using the normal ray projection method, and the CT value H(p) corresponding to the point is extracted. The thickness and CT values ​​are normalized to obtain the normalized thickness. and normalized CT value ; The stiffness exponent S(p) of the point is calculated using the following formula: in, and These are the preset weighting coefficients.

3. The system according to claim 1, characterized in that, The method of combining vascular diameter and lung lobe functional density distribution to assign risk value to spatial voxels and construct a surgical risk potential field includes: For any voxel v in space, the risk cost R(v) of the voxel is determined by the following formula: in, The stiffness index is the stiffness index of the anatomical barrier closest to voxel v. The diameter of the blood vessel closest to voxel v. This represents the functional density of the lung lobe at the location of voxel v. , These are the dimensionless values ​​of each term after min-max normalization. , , These are positive weighting coefficients.

4. The system according to claim 1, characterized in that, The calculation of the surface irregularity of the lesion determines the multi-point target region on the lesion surface, including: Calculate the average curvature of each vertex on the surface mesh of the 3D model of the lesion; The set of vertices with an average curvature greater than a preset threshold is defined as candidate target points; The candidate target points are clustered into 5 clusters using the K-means clustering algorithm, and the centroid of each cluster is used as a target point in the multi-point target region.

5. The system according to claim 1, characterized in that, The iteration step size and direction of the search unit are controlled by the local gradient tensor of the risk potential field, including: At the current position of the search unit, calculate the gradient tensor G of the risk potential field at that point; Perform eigenvalue decomposition on the gradient tensor G to obtain eigenvalues. and the corresponding feature vector ; will be with the minimum eigenvalue Corresponding feature vector The priority iteration direction for the search unit; the iteration step size L is set by the following formula: in, The base step size is k, and the adjustment coefficient is k.

6. The system according to claim 1, characterized in that, The process of generating a variable-boundary safe passage around the initial path based on the risk potential gradient around the initial path and the local curvature includes: For each point p on the initial path, calculate the local curvature at that point. and the mean gradient of the risk potential field in the local region surrounding the point ; The radial width W(p) of the channel at the point is calculated using the following formula: in, Based on the width, and It is a positive proportionality coefficient; A circular cross section is constructed on the normal plane with the radial width W(p) of each point on the path, and all the circular cross sections are connected to form the safety passage.

7. The system according to claim 1, characterized in that, The step of smoothing and optimizing the initial path using the boundary of the safety passage as a constraint to obtain the planned path includes: The initial path is represented as a series of discrete path points. ; Construct the target energy function: Among them, the smoothing term The integral of the total curvature of the path, and the boundary constraint terms. The penalty function is the distance from the path point to the safety passage boundary, and the value of the function increases when the path point exceeds the passage boundary; The gradient descent method is used to iteratively adjust the positions of the path points until the target energy function E(P) converges to the minimum value, thus obtaining the planned path.

Images