3D Image Reconstruction and Localization System and Method for Lung Puncture

Abstract

The present invention relates to the technical field of image reconstruction, and particularly to a three-dimensional imaging reconstruction positioning system and method for lung puncture. The method includes the following steps: obtaining CT scan sequence data of a patient and multi-point infrared marker data of the chest wall, where the multi-point infrared marker data of the chest wall includes respiratory reference position data, real-time angular velocity data of the chest wall motion state, and real-time acceleration data of each point of the chest wall; based on the multi-point infrared marker data of the chest wall, a non-linear viscoelastic prediction model of lung tissue based on the Maxwell model is established, and the real-time angular velocity data and the real-time acceleration data are input into the non-linear viscoelastic prediction model of lung tissue for iterative calculation to obtain lung tissue stress-strain relationship data. The present invention can simplify the non-feature region while maintaining the grid density of the feature region, greatly reducing the computational burden and ensuring the details and global consistency of the reconstruction.

Description

Technical Field

[0001] The present invention relates to the technical field of image reconstruction, and in particular to a three-dimensional image reconstruction and positioning system and method for lung puncture. Background Art

[0002] Lung puncture refers to an operation of obtaining lung tissue or fluid samples through a puncture needle for diagnosis or treatment; percutaneous lung puncture is a specific way of lung puncture, specifically referring to performing biopsy or other operations by puncturing through the skin, passing through the chest wall and pleura into the lung tissue under the guidance of imaging equipment (such as CT, B-ultrasound, X-ray, etc.). Three-dimensional image reconstruction is to process the CT image data (usually in DICOM format) of the patient's lungs and use professional software to reconstruct the three-dimensional anatomical structure of the lungs. Through accurate three-dimensional positioning, the one-time puncture success rate is significantly improved. For example, when using a robotic three-dimensional targeted positioning system for CT-guided lung biopsy, the one-time puncture success rate can reach 75%, which is much higher than the traditional method.

[0003] However, traditional three-dimensional image reconstruction and positioning methods for lung puncture often have the following problems: During lung tumor biopsy, doctors need to accurately locate the tumor position and obtain tissue samples. Three-dimensional reconstruction technology can help operators plan the best puncture path and avoid important blood vessels and tracheas. However, the technical difficulty in this case lies in that respiratory movement will cause the tumor position to change, making real-time tracking difficult; if the tumor is small (<1 cm) or located deep, the positioning accuracy requirement is extremely high, increasing the technical difficulty. Summary of the Invention

[0004] Based on this, it is necessary for the present invention to provide a three-dimensional image reconstruction and positioning system and method for lung puncture to solve at least one of the above technical problems.

[0005] To achieve the above object, a three-dimensional image reconstruction and positioning method for lung puncture includes the following steps:

[0006] Step S1: Obtain the CT scan sequence data of the patient and the multi-point infrared marker data of the chest wall, where the multi-point infrared marker data of the chest wall includes respiratory reference position data, real-time angular velocity data of the chest wall movement state, and real-time acceleration data of each point on the chest wall;

[0007] Step S2: Based on the multi-point infrared marker data of the chest wall, a non-linear viscoelastic prediction model of lung tissue based on the Maxwell model, and input the real-time angular velocity data and real-time acceleration data into the non-linear viscoelastic prediction model of lung tissue for iterative calculation to obtain lung tissue stress-strain relationship data, where the stress-strain relationship data includes instantaneous stress response data and delayed strain component data;

[0008] Step S3: Based on the stress-strain relationship data, perform partitioned mesh generation on the lung tissue to obtain the first deformation field data of the central region based on high-density meshes, the second deformation field data of the peripheral region based on low-density meshes, and the third deformation field data of the transition region constructed with adaptive meshes based on local stress gradients in the remaining region;

[0009] Step S4: Use the optical flow method to calculate the motion vector fields of the first deformation field data of the central region, the second deformation field data of the peripheral region, and the third deformation field data of the transition region, construct an energy functional according to the brightness constancy constraint and the velocity smoothness constraint, and solve it by the variational method to obtain the optimized tissue deformation field data;

[0010] Step S5: Apply the tissue deformation field data to the CT scan sequence data and perform 3D reconstruction based on the bubble animation simulation algorithm, simplify the meshes in the non-feature line regions while maintaining the mesh density in the feature line regions, to obtain the puncture path planning model data with respiratory phase mapping.

[0011] The present invention provides comprehensive and accurate imaging and kinematic information for subsequent analysis by acquiring the CT scan sequence data and multi-point infrared marker data of the chest wall of the patient. Among them, the introduction of the respiratory reference position, real-time angular velocity, and acceleration data enables the precise capture of the dynamic changes of the lung tissue during the respiratory cycle, laying a foundation for subsequent modeling and deformation analysis. The fusion of such multi-dimensional data effectively improves the accuracy and reliability of lung puncture positioning and avoids puncture deviation caused by respiratory movement. The construction of a non-linear viscoelastic prediction model of lung tissue based on the Maxwell model is one of the core innovations of this method. By inputting the real-time angular velocity and acceleration data into the model for iterative calculation, the stress-strain relationship of the lung tissue during respiration can be accurately predicted. This prediction method based on a physical model not only considers the non-linear viscoelastic characteristics of the lung tissue but also optimizes the accuracy of the model through the real-time feedback of dynamic data. The acquisition of the instantaneous stress response data and the delayed strain component data provides key mechanical parameters for subsequent mesh generation and deformation field calculation, enabling the puncture path planning to better conform to the actual physiological state of the lung tissue. The zonal mesh generation further optimizes the computational efficiency and accuracy. By dividing the lung tissue into a central region, a peripheral region, and a transition region, and respectively using high-density meshes, low-density meshes, and adaptive meshes, the deformation characteristics of different regions can be accurately captured. The high-density meshes in the central region ensure accurate modeling near the puncture target point, the low-density meshes in the peripheral region reduce the computational burden, and the adaptive meshes in the transition region can flexibly respond to regions with large stress gradients. This zonal strategy not only improves the flexibility of mesh generation but also effectively balances the allocation of computational resources, making the entire modeling process more efficient. The application of the optical flow method further improves the accuracy of tissue deformation field calculation. By constructing an energy functional by combining the brightness constancy constraint and the velocity smoothness constraint and solving it using the variational method, the noise and discontinuity in the motion vector field can be effectively eliminated. The optimized tissue deformation field data can not only accurately reflect the dynamic changes of the lung tissue during respiration but also provide high-quality input data for subsequent 3D reconstruction. This motion vector field calculation method based on the optical flow method significantly improves the computational accuracy and stability of the deformation field and provides a reliable basis for puncture path planning. By applying the tissue deformation field data to the CT scan sequence data and combining it with the bubble animation simulation algorithm for 3D reconstruction, a puncture path planning model with respiratory phase mapping is constructed. During the reconstruction process, while maintaining the mesh density in the feature line region, the mesh in the non-feature line region is simplified, which not only retains the detailed information of important structures but also further optimizes the complexity of the model. This 3D reconstruction method can intuitively display the dynamic changes of the lung tissue at different respiratory phases, providing clear puncture path planning and real-time navigation support for clinicians and effectively improving the success rate and safety of puncture.In summary, through multi-dimensional data acquisition, precise mechanical modeling, efficient mesh generation, accurate calculation of the motion vector field, and optimized 3D reconstruction, this method significantly improves the positioning accuracy and safety of lung puncture, provides strong technical support for clinical lung puncture surgery, and has broad application prospects and important clinical significance.

[0012] Preferably, the present invention further provides a 3D image reconstruction and positioning system for lung puncture, which is used to execute the above-mentioned 3D image reconstruction and positioning method for lung puncture. The 3D image reconstruction and positioning system for lung puncture includes:

[0013] A data acquisition module, which is used to obtain the CT scan sequence data of the patient and the multi-point infrared marker data of the chest wall. The multi-point infrared marker data of the chest wall includes the respiratory reference position data, the real-time angular velocity data of the chest wall motion state, and the real-time acceleration data of each point on the chest wall;

[0014] A tissue mechanics modeling module, which is used to predict the non-linear viscoelasticity of lung tissue based on the multi-point infrared marker data of the chest wall using the Maxwell model, and input the real-time angular velocity data and the real-time acceleration data into the non-linear viscoelasticity prediction model of lung tissue for iterative calculation to obtain the stress-strain relationship data of lung tissue. The stress-strain relationship data includes the instantaneous stress response data and the delayed strain component data;

[0015] A partition mesh generation module, which is used to perform partition mesh generation on the lung tissue based on the stress-strain relationship data to obtain the first deformation field data of the central region based on high-density meshes, the second deformation field data of the peripheral region based on low-density meshes, and the third deformation field data of the transition region constructed by using adaptive meshes based on local stress gradients in the remaining region;

[0016] A motion vector calculation module, which is used to calculate the motion vector field of the first deformation field data of the central region, the second deformation field data of the peripheral region, and the third deformation field data of the transition region by using the optical flow method, construct an energy functional according to the brightness constancy constraint and the velocity smoothness constraint, and solve it by variational method to obtain the optimized tissue deformation field data;

[0017] A 3D reconstruction and path planning module, which is used to apply the tissue deformation field data to the CT scan sequence data and perform 3D reconstruction based on the bubble animation simulation algorithm, simplify the meshes in the non-feature line region while maintaining the mesh density in the feature line region, and obtain the puncture path planning model data with respiratory phase mapping.

[0018] In the present invention, the data acquisition module provides high-quality basic information for subsequent modeling and analysis by acquiring the CT scan sequence data of the patient and the multi-point infrared marker data of the chest wall. The CT scan sequence data can record the anatomical structure of the lungs in detail, while the multi-point infrared marker data of the chest wall accurately captures the impact of respiratory motion on lung tissue through the respiratory reference position, real-time angular velocity, and acceleration information. The fusion of such multi-dimensional data enables the system to comprehensively consider the dynamic changes of lung tissue during the respiratory cycle, providing accurate real-time data support for puncture path planning. The tissue mechanics modeling module constructs a non-linear viscoelastic prediction model of lung tissue based on the Maxwell model and performs iterative calculations through real-time angular velocity and acceleration data. This process not only considers the complex mechanical properties of lung tissue but also optimizes the accuracy of the model through the real-time feedback of dynamic data. The generation of instantaneous stress response and delayed strain component data provides key mechanical parameters for subsequent deformation field calculations, enabling the system to accurately simulate the dynamic deformation behavior of lung tissue during respiration and providing a reliable mechanical basis for puncture path planning. The partitioned mesh generation module performs partitioned mesh generation on lung tissue according to stress-strain relationship data, generating deformation field data for the central region, peripheral region, and transition region respectively. This partitioning strategy achieves a balance between computational accuracy and efficiency by focusing on the key puncture regions with high-density meshes, covering the peripheral regions with low-density meshes, and optimizing the transition regions with adaptive meshes based on local stress gradients. Partitioned mesh generation not only improves the accuracy of deformation field calculations but also provides an efficient mesh basis for subsequent 3D reconstruction and path planning, ensuring the scientificity and feasibility of puncture path planning. The motion vector calculation module uses the optical flow method combined with the brightness constancy constraint and the velocity smoothness constraint to calculate the motion vector field for the deformation field data of each region. By constructing an energy functional and solving it using the variational method, the system can optimize the motion vector field and eliminate noise and discontinuities. The optimized tissue deformation field data can accurately reflect the dynamic changes of lung tissue during respiration, providing high-quality dynamic information support for puncture path planning and further enhancing the accuracy and safety of the puncture path. Finally, the 3D reconstruction and path planning module applies the tissue deformation field data to the CT scan sequence data and performs 3D reconstruction through the bubble animation simulation algorithm. During the reconstruction process, the system maintains the mesh density in the feature line regions while simplifying the meshes in the non-feature line regions, generating a puncture path planning model with respiratory phase mapping. This 3D reconstruction strategy not only retains the detailed information of important anatomical structures but also ensures that the puncture path can adapt to the dynamic changes of lung tissue through respiratory phase mapping. The finally generated puncture path planning model can provide intuitive and reliable navigation support for clinicians, significantly improving the success rate and safety of puncture surgery. Through the collaborative work of the above modules, the lung puncture 3D image reconstruction and positioning system realizes the full-process optimization from data acquisition to puncture path planning.Each module not only complements each other technically, but also forms an organic whole functionally, providing precise and efficient technical support for lung puncture surgery and promoting the development of clinical puncture technology to a higher level. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Other features, objects, and advantages of the present invention will become more apparent from the following detailed description of non-limiting embodiments read with reference to the accompanying drawings:

[0020] Figure 1 It is a schematic flow chart of the steps of the three-dimensional image reconstruction and positioning method for lung puncture of the present invention;

[0021] Figure 2 is Figure 1 a detailed schematic flow chart of step S1 in

[0022] Figure 3 is Figure 1 a detailed schematic flow chart of step S2 in DETAILED DESCRIPTION OF THE EMBODIMENTS

[0023] The technical method of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of them. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative efforts belong to the scope of protection of the present invention.

[0024] In addition, the drawings are only schematic diagrams of the present invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and thus their repeated description will be omitted. Some of the block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. The functional entities can be implemented in software form, or in one or more hardware modules or integrated circuits, or in different networks and / or processor methods and / or microcontroller methods.

[0025] It should be understood that although the terms "first", "second", etc. may be used here to describe each unit, these units should not be limited by these terms. These terms are only used to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, the first unit can be called the second unit, and similarly the second unit can be called the first unit. The term "and / or" used here includes any and all combinations of one or more of the listed associated items.

[0026] To achieve the above object, please refer to Figures 1 to 3, the present invention provides a three-dimensional image reconstruction and positioning method for lung puncture, and the method includes the following steps:

[0027] Step S1: Obtain the CT scan sequence data of the patient and the multi-point infrared marker data of the chest wall, where the multi-point infrared marker data of the chest wall includes the respiratory reference position data, the real-time angular velocity data of the chest wall movement state, and the real-time acceleration data of each point on the chest wall;

[0028] In an embodiment of the present invention, a high-resolution spiral CT scanner is used to perform multi-layer scans on the patient's chest to collect the CT image sequence data of the patient in different respiratory states. The image resolution is 0.5mm×0.5mm, the slice interval is 1mm, covering the entire lung range. At the same time, high-precision infrared marker points are arranged in multiple areas of the patient's chest wall. Each marker point collects position information, angular velocity data, and acceleration data in real time through an infrared tracker. The real-time acceleration data refers to the acceleration change of each point on the chest wall. Since the movement of the chest wall is closely related to the expansion and contraction of the lungs, this acceleration data represents the linear acceleration experienced by each point on the chest wall surface during breathing; the real-time angular velocity data refers to the angular velocity of the chest wall movement during breathing (i.e., the rotation rate of the chest wall in three-dimensional space) as the lungs expand and contract. This data can be obtained by tracking the rotation of the chest wall with infrared markers, reflecting the degree of rotation or swing of each point on the chest wall. The acquisition frequency is set to 120Hz to ensure that the dynamic changes in the rapid breathing state can be captured. The infrared marker point data includes the reference position data for calibrating the initial respiratory state, and the angular velocity and acceleration data for capturing the dynamic respiratory deformation information. These data are synchronized through a dedicated data acquisition module and stored in a standardized format for subsequent processing.

[0029] Step S2: Based on the multi-point infrared marker data of the chest wall, a non-linear viscoelastic prediction model of lung tissue based on the Maxwell model is established, and the real-time angular velocity data and the real-time acceleration data are input into the non-linear viscoelastic prediction model of lung tissue for iterative calculation to obtain the stress-strain relationship data of lung tissue, where the stress-strain relationship data includes the instantaneous stress response data and the delayed strain component data;

[0030] In an embodiment of the present invention, the non-linear viscoelastic characteristics of the Maxwell model are used to describe the mechanical behavior of lung tissue. First, based on the collected infrared marker data, a prediction model of lung tissue based on dynamic stress-strain behavior is constructed. The model formula is: ; where is the elastic modulus, is the viscous modulus, is the viscous delay time, and the initial value is obtained through experimental calibration . Taking the real-time angular velocity data and acceleration data as inputs, combining with the patient's reference respiratory cycle, and using the nonlinear finite element method for iterative calculation with each iteration step set to 0.01 s, the instantaneous stress response data and the delayed strain component data are obtained, forming a time-series stress-strain relationship of the lung tissue.

[0031] Step S3: Based on the stress-strain relationship data, perform partitioned mesh generation for the lung tissue to obtain the first deformation field data of the central region based on high-density meshes, the second deformation field data of the peripheral region based on low-density meshes, and the third deformation field data of the transition region constructed by using adaptive meshes based on local stress gradients in the remaining region;

[0032] In the embodiment of the present invention, when using the stress-strain relationship data for lung tissue mesh generation, the entire lung is first divided into a central region, a peripheral region, and a transition region according to the anatomical characteristics and stress distribution of the lung tissue. The central region adopts high-density regular meshes, and the size of each mesh unit is set to 0.2 mm × 0.2 mm to ensure the deformation calculation accuracy of this region; the peripheral region adopts low-density regular meshes, and the mesh unit size is set to 0.5 mm × 0.5 mm; the transition region generates adaptive meshes through local stress gradients, and the mesh density is dynamically adjusted according to the stress gradient, with the minimum mesh unit size being 0.3 mm × 0.3 mm and the highest density being the same as that of the central region. The mesh generation adopts a voxelization method, realizes three-dimensional reconstruction by layer-by-layer scanning of CT data, and updates the dynamic deformation in combination with the patient's respiratory state in real time.

[0033] Step S4: Use the optical flow method to calculate the motion vector fields for the first deformation field data of the central region, the second deformation field data of the peripheral region, and the third deformation field data of the transition region, and construct an energy functional according to the brightness constancy constraint and the velocity smoothness constraint, and solve through the variational method to obtain the optimized tissue deformation field data;

[0034] In the embodiment of the present invention, the optical flow method is used to calculate the motion vector fields for the deformation field data in the partitioned meshes, and the basic optical flow constraint formula is used, where and are the motion vector components. Calculate the first deformation field data of the central region, the second deformation field data of the peripheral region, and the third deformation field data of the transition region separately, and construct an energy functional in combination with the brightness constancy constraint and the velocity smoothness constraint: . Solve this energy functional through the variational method, and use the gradient descent optimization algorithm to iteratively solve to obtain the optimized tissue deformation field data. The parameter in the brightness constancy constraint is set to to balance the influence of brightness and velocity.

[0035] Step S5: Apply the tissue deformation field data to the CT scan sequence data and perform three-dimensional reconstruction based on the bubble animation simulation algorithm. While maintaining the grid density in the feature line region, simplify the grid in the non-feature line region to obtain the puncture path planning model data with respiratory phase mapping.

[0036] In the embodiment of the present invention, the optimized tissue deformation field data is applied to the CT scan sequence data, and three-dimensional reconstruction of lung deformation is performed through the bubble animation simulation algorithm. The bubble animation algorithm is driven by the dynamic data of the patient's respiratory phase, maintains a high-density grid (cell size 0.2mm×0.2mm×0.2mm) in the feature line region, and reduces the computational complexity in the non-feature line region through grid simplification. The simplified grid density is 40% of the original density. During the three-dimensional reconstruction process, the respiratory phase mapping technology is used to generate corresponding respiratory state labels for each reconstructed three-dimensional shape, and finally, three-dimensional model data available for puncture path planning is generated, and the path planning accuracy is ensured to reach the sub-millimeter level through simulation verification.

[0037] The present invention provides comprehensive and accurate imaging and kinematic information for subsequent analysis by acquiring the CT scan sequence data and multi-point infrared marker data of the chest wall of the patient. Among them, the introduction of the respiratory reference position, real-time angular velocity and acceleration data enables the precise capture of the dynamic changes of the lung tissue during the respiratory cycle, laying a foundation for subsequent modeling and deformation analysis. The fusion of such multi-dimensional data effectively improves the accuracy and reliability of lung puncture positioning, and avoids puncture deviation caused by respiratory movement. The construction of a non-linear viscoelastic prediction model of lung tissue based on the Maxwell model is one of the core innovations of this method. By inputting the real-time angular velocity and acceleration data into the model for iterative calculation, the stress-strain relationship of the lung tissue during respiration can be accurately predicted. This prediction method based on a physical model not only considers the non-linear viscoelastic characteristics of the lung tissue, but also optimizes the accuracy of the model through the real-time feedback of dynamic data. The acquisition of the instantaneous stress response data and the delayed strain component data provides key mechanical parameters for subsequent mesh generation and deformation field calculation, enabling the puncture path planning to be more in line with the actual physiological state of the lung tissue. The zonal mesh generation further optimizes the computational efficiency and accuracy. By dividing the lung tissue into a central region, a peripheral region and a transition region, and respectively using high-density meshes, low-density meshes and adaptive meshes, the deformation characteristics of different regions can be accurately captured. The high-density meshes in the central region ensure accurate modeling near the puncture target point, the low-density meshes in the peripheral region reduce the computational burden, and the adaptive meshes in the transition region can flexibly respond to regions with large stress gradients. This zonal strategy not only improves the flexibility of mesh generation, but also effectively balances the allocation of computational resources, making the entire modeling process more efficient. The application of the optical flow method further improves the accuracy of tissue deformation field calculation. By constructing an energy functional by combining the brightness constancy constraint and the velocity smoothness constraint, and solving it using the variational method, the noise and discontinuity in the motion vector field can be effectively eliminated. The optimized tissue deformation field data can not only accurately reflect the dynamic changes of the lung tissue during respiration, but also provide high-quality input data for subsequent 3D reconstruction. This motion vector field calculation method based on the optical flow method significantly improves the calculation accuracy and stability of the deformation field, providing a reliable basis for puncture path planning. By applying the tissue deformation field data to the CT scan sequence data and combining with the bubble animation simulation algorithm for 3D reconstruction, a puncture path planning model with respiratory phase mapping is constructed. During the reconstruction process, while maintaining the mesh density in the feature line region, the mesh in the non-feature line region is simplified, which not only retains the detailed information of important structures, but also further optimizes the complexity of the model. This 3D reconstruction method can intuitively display the dynamic changes of the lung tissue at different respiratory phases, providing clear puncture path planning and real-time navigation support for clinicians, and effectively improving the success rate and safety of puncture.In summary, through multi-dimensional data acquisition, precise mechanical modeling, efficient mesh generation, accurate calculation of the motion vector field, and optimized 3D reconstruction, this method significantly improves the positioning accuracy and safety of lung puncture, provides strong technical support for clinical lung puncture surgery, and has broad application prospects and important clinical significance.

[0038] Preferably, step S1 includes the following steps:

[0039] Step S11: Acquire CT scan sequence data during the patient's respiratory cycle;

[0040] In the embodiment of the present invention, the patient lies in the supine position, and is scanned using a high-resolution spiral CT scanner. The scanning range covers the upper edge of the chest cavity to the lower edge of the abdominal cavity to ensure that the complete lung area is included. The CT scan is performed with a spatial resolution of 0.5mm×0.5mm, and the thickness of each scanned layer is 1mm. During the acquisition process, synchronous electrocardiogram-respiratory gating technology is adopted, and the patient's respiratory cycle is recorded by an external monitoring device and the respiratory phase state (such as end-inspiration, end-expiration, etc.) of the CT image sequence is synchronously marked. The exposure parameters of the CT device are 120 kVp and the current is 200 mAs to ensure image quality. At the same time, to avoid motion artifacts, the patient is required to hold their breath for a short time at a specific respiratory phase.

[0041] Step S12: Arrange a plurality of infrared reflection marker points on the surface of the patient's chest wall to obtain infrared reflection marker point data, where a marker point is arranged at the manubrium sterni position to obtain respiratory reference position information, marker points are arranged at the ends of the left and right seventh ribs respectively to obtain angular velocity information, and a marker point is arranged 5 cm above the umbilicus of the abdomen to obtain acceleration information;

[0042] In the embodiment of the present invention, the respiratory dynamic information of the patient is obtained in real time through infrared reflection marker points. First, spherical reflection marker points with a diameter of 5 mm are selected, and their surfaces are coated with a high-reflectivity coating to ensure the accuracy of optical capture. A marker point is arranged at the manubrium sterni position to record the respiratory reference position information of the patient; marker points are arranged at the ends of the left and right seventh ribs respectively to capture the dynamic angular velocity changes of the chest wall during the patient's breathing; a marker point is arranged 5 cm above the umbilicus of the abdomen to measure the longitudinal acceleration caused by breathing. The arrangement of the marker points is fixed by a special medical-grade tape to ensure that the stability is not disturbed by skin movement during breathing.

[0043] Step S13: Dynamically track the infrared reflection marker point data using a high-speed infrared camera system, and obtain the three-dimensional spatial coordinate data of each marker point through an image registration algorithm;

[0044] In the embodiment of the present invention, a high-speed infrared imaging system (sampling frequency: 120 Hz, resolution: 1920×1080 pixels) is adopted to dynamically track the infrared reflection marker points arranged on the patient's chest wall. In actual operation, the imaging system is installed on a bracket 1 meter above the patient, and the camera angle is adjusted to cover all marker point areas. The system irradiates the marker points with a laser through an infrared light source, and at the same time uses a high-sensitivity infrared sensor to capture the reflection signal. Through an image registration algorithm (such as a registration algorithm based on SURF feature point matching), the two-dimensional pixel coordinates of each marker point are extracted from consecutive frame images, and the three-dimensional spatial coordinates of the marker points are inversely calculated through the calibration parameters (focal length, distortion coefficient, etc.) of the camera.

[0045] Step S14: Smooth the motion trajectory of the marker points based on the three-dimensional spatial coordinate data, so as to obtain the respiration reference position data with the influence of noise eliminated.

[0046] In the embodiment of the present invention, by analyzing the collected three-dimensional coordinate data of the infrared marker points, it is found that the original trajectory contains high-frequency noise caused by equipment errors and environmental interference. To eliminate the noise, a trajectory smoothing algorithm based on quintic spline interpolation is adopted to perform curve fitting on the original data with a sampling frequency of 120 Hz. The specific steps of the smoothing process include: 1) fitting independent spline curves for the coordinates of each marker point respectively, 2) selecting an appropriate smoothing coefficient (such as 0.05) to ensure that both the noise can be eliminated and the key motion information is not lost, 3) resampling the fitted trajectory curve point by point to generate the smoothed respiration reference position data.

[0047] Step S15: Perform a derivative operation on the motion trajectory of the marker points according to the three-dimensional spatial coordinate data to obtain the real-time angular velocity data and the real-time acceleration data respectively.

[0048] In the embodiment of the present invention, the smoothed three-dimensional trajectory data is used to calculate the motion parameters of the marker points by using the numerical differentiation method. The specific operation is as follows: the first-order derivative of the trajectory data is calculated by the central difference method to obtain the real-time angular velocity data of each marker point; then the acceleration data is obtained through the second-order derivative calculation. Taking the marker points at the ends of the seventh ribs on the left and right as an example, the real-time angular velocity is obtained by calculating the angular displacement change at each moment and combining the sampling frequency; for the abdominal marker points, the acceleration is obtained through the second-order derivative of the longitudinal displacement. To verify the accuracy of the derivative results, the calculated angular velocity and acceleration are compared with the theoretical respiration cycle model, and the error is controlled within 5%.

[0049] Step S16: Align the respiration reference position data, the real-time angular velocity data, and the real-time acceleration data in time series, so as to obtain the multi-point infrared marker data of the chest wall.

[0050] In the embodiments of the present invention, since there may be a deviation in the time axes of CT image acquisition and infrared marker data acquisition, timing alignment processing needs to be performed. The specific method is as follows: First, extract the respiratory reference position curve in the infrared marker data, perform similarity matching with the respiratory gating signal of CT scanning, and use an algorithm based on dynamic time warping (DTW) to align the two. After alignment, according to the unified time axis, combine the respiratory reference position data, real-time angular velocity data, and real-time acceleration data into multi-point chest wall infrared marker data, and store it in a standardized data format indexed by timestamp (such as CSV or HDF5 file). The finally aligned multi-point infrared marker data is used for subsequent lung tissue deformation modeling and analysis.

[0051] The present invention collects CT scan sequence data within the patient's respiratory cycle, which can completely record the morphological changes of the lungs at different respiratory stages. This dynamic data collection method provides a rich imaging basis for subsequent analysis, ensuring that the puncture path planning can fully consider the influence of respiratory movement, thereby improving the accuracy and safety of puncture. By arranging multiple infrared reflection marker points on the patient's chest wall surface, multi-dimensional monitoring of respiratory movement is achieved. Marker points are respectively arranged at the manubrium sterni, the end of the seventh rib, and above the umbilicus of the abdomen, enabling accurate acquisition of respiratory reference position, angular velocity, and acceleration information. This layout makes full use of human anatomical features to ensure the comprehensiveness and accuracy of data collection. A high-speed infrared camera system is used to dynamically track the infrared reflection marker points, and the three-dimensional spatial coordinate data of each marker point is obtained through an image registration algorithm. The high-speed camera system can capture the movement trajectory of the marker points in real time, while the image registration algorithm further improves the accuracy and reliability of the data, providing high-precision three-dimensional motion data for subsequent analysis. By smoothing the movement trajectory of the marker points, the influence of noise on the respiratory reference position data is eliminated. This data preprocessing method effectively improves the signal-to-noise ratio of the signal, making the respiratory reference position data more stable and reliable, and laying a foundation for subsequent angular velocity and acceleration calculations. Real-time angular velocity and acceleration data are extracted from the three-dimensional spatial coordinate data through derivative operations. These data reflect the dynamic characteristics of respiratory movement and can provide key kinematic parameters for the deformation analysis of lung tissue. This derivative calculation method based on time series ensures the real-time and dynamic nature of the data, enabling subsequent models to accurately simulate the movement state of lung tissue. By performing temporal alignment on the respiratory reference position data, real-time angular velocity data, and real-time acceleration data, a complete multi-point infrared marker data of the chest wall is integrated. This data alignment method ensures the synchronization between different data types, enabling subsequent models to perform analysis and calculations based on a unified time reference, further improving the overall quality and usability of the data. In summary, through a refined data collection and processing process, this method can efficiently and accurately obtain the dynamic information of the lungs within the patient's respiratory cycle. These high-quality data provide a solid foundation for subsequent lung tissue deformation analysis and puncture path planning, significantly improving the accuracy and safety of lung puncture surgery, and having important clinical application value.

[0052] Preferably, step S2 includes the following steps:

[0053] Step S21: Obtain the elastic modulus parameter data and viscosity coefficient parameter data of the standard human lung tissue;

[0054] In the embodiments of the present invention, by referring to the literature and experimental data, the elastic modulus and viscous coefficient parameter data of the standard human lung tissue are obtained. Specifically, in the laboratory, the excised lung tissue samples are used, and their elastic modulus (unit: kPa) is obtained through a tensile test, and their viscous coefficient (unit: Pa·s) is measured through a shear test. For example, a lung tissue sample from the lower lobe of the left lung of a healthy adult is selected and tested in an environment with a room temperature of 22°C and a humidity of 60%. The measured result of the elastic modulus is 4.2±0.3 kPa, and the measured result of the viscous coefficient is 2.1±0.2 Pa·s. In addition, to further verify the data, a standard simulation software (such as ANSYS or COMSOL) is used to verify the above parameters through a model to ensure its applicability to the subsequent modeling process.

[0055] Step S22: Establish an elastic element model in the Maxwell model according to the elastic modulus parameter data of the lung tissue, so as to obtain initial elastic response data; establish a viscous element model in the Maxwell model according to the viscous coefficient parameter data of the lung tissue, so as to obtain initial viscous response data;

[0056] In the embodiments of the present invention, the obtained elastic modulus and viscous coefficient parameters are used to establish an elastic element model and a viscous element model in the Maxwell model respectively. The elastic element model is based on Hooke's law, and the initial elastic response data is obtained by calculating the product of the elastic modulus and the strain. The viscous element model is based on Newton's viscosity law, and the initial viscous response data is obtained by calculating the product of the viscous coefficient and the strain rate. In the specific implementation of the calculation, the MATLAB tool is used to input the elastic modulus and viscous coefficient data respectively, and the simulation calculation is carried out in combination with a discrete time step of 0.01 s. For example, under the condition of applying a constant strain of 0.5%, the calculated elastic response is 0.021 kPa, and the viscous response is 0.0105 Pa·s.

[0057] Step S23: Combine the initial elastic response data and the initial viscous response data in series to obtain the basic Maxwell unit model data;

[0058] In the embodiments of the present invention, the elastic response data and the viscous response data obtained in step S22 are combined in series to establish a basic Maxwell unit model. In the specific implementation, combined with the series relationship of the Maxwell model, that is, the stress is equal to the sum of the stresses of the elastic and viscous elements, and the strain rate is equal to the strain rate of the viscous element minus the strain rate of the elastic element, the dynamic coupling response data of the two are solved by a numerical method. A simulation program is written in Python, and the data matrices of the elastic response and the viscous response are input, and the curve of the stress changing with time is obtained through calculation. For example, within 1 second, the stress response curve of the basic Maxwell unit model obtained by calculation shows the characteristics of rising rapidly and then gradually stabilizing.

[0059] Step S24: Establish a multi-dimensional state equation based on the basic Maxwell element model data and the multi-point infrared marker data of the chest wall, and perform numerical discretization to obtain a non-linear viscoelastic prediction model of lung tissue;

[0060] In the embodiment of the present invention, the basic Maxwell element model data is combined with the multi-point infrared marker data of the chest wall (including real-time angular velocity and acceleration data) to construct a multi-dimensional state equation. The state equation is based on the non-linear viscoelastic theory and has the form of: ; where is the stress, is the elastic modulus, is the viscosity coefficient, is the strain. By introducing the angular velocity and acceleration data of the infrared marker points, the time-related terms of the state equation are established. During the calculation process, the finite difference method is used to discretize the continuous equation. In a specific implementation, based on a time step of 0.005 seconds, a simulation environment is built using MATLAB Simulink, and the basic Maxwell element model data and the multi-point infrared marker data of the chest wall are input to obtain a non-linear viscoelastic prediction model of lung tissue. In the model test, with a strain rate of 0.2 / s for the assumed initial conditions, the deviation between the simulated state equation response and the real data is controlled within 3%.

[0061] Step S25: Input the real-time angular velocity data and the real-time acceleration data into the non-linear viscoelastic prediction model of lung tissue for iterative calculation to obtain the stress-strain relationship data of lung tissue, where the stress-strain relationship data includes instantaneous stress response data and delayed strain component data.

[0062] In the embodiment of the present invention, the non-linear viscoelastic prediction model of lung tissue constructed in Step S24 is used, and the real-time angular velocity data and acceleration data are used as inputs to obtain the stress-strain relationship data through iterative solution. The specific operations are as follows: 1) Taking the angular velocity and acceleration data of the chest wall infrared marker points as inputs, calculate the instantaneous stress response; 2) Calculate the delayed strain component according to the delayed characteristics of the viscous element. Use Python and the NumPy library to write code, set the number of iterative steps to 200 times, and the initial conditions are that both the stress and strain are zero. Run to obtain a stress-strain data matrix with a time resolution of 0.01 seconds. The calculation results show that the instantaneous stress response data has the characteristics of high-frequency fluctuations, and the delayed strain component data shows exponential decay over time. This stress-strain relationship data is finally used for mesh generation and deformation field calculation in the lung tissue region.

[0063] The present invention obtains the elastic modulus and viscosity coefficient parameter data of standard human lung tissue, and provides key physical property parameters for the establishment of subsequent models. These parameters reflect the mechanical behavior of lung tissue under physiological conditions and are the basis for achieving accurate modeling. By introducing standard parameter data, the universality and accuracy of the model can be ensured, making it suitable for the analysis of lung tissue characteristics of different patients. Based on the elastic modulus and viscosity coefficient parameters of lung tissue, the elastic element and viscosity element models in the Maxwell model are respectively established, and the initial elastic response and viscosity response data are obtained. This step-by-step modeling method can clearly simulate the elastic deformation and viscous flow characteristics of lung tissue when subjected to force, and provides an accurate mechanical response basis for subsequent comprehensive modeling. By modeling the elastic and viscous characteristics of lung tissue separately, its complex mechanical behavior can be more comprehensively reflected. The initial elastic response data and the viscosity response data are combined in series to form the basic Maxwell unit model data. This series connection method simulates the nonlinear viscoelastic behavior of lung tissue under dynamic loading, so that the model can be more in line with the actual physiological state. The establishment of the basic Maxwell unit model provides a core structure for the subsequent multidimensional state equation, ensuring the applicability of the model in a complex mechanical environment. Combining the basic Maxwell unit model data and the chest wall multi-point infrared marker data, a multi-dimensional state equation was established, and the nonlinear viscoelastic prediction model of lung tissue was obtained through numerical discretization. This process combines the mechanical properties of lung tissue with the actual respiratory motion data, so that the model can dynamically predict the stress-strain relationship of lung tissue during breathing. Through numerical discretization, the computational efficiency and adaptability of the model are further improved, so that it can reflect the dynamic changes of lung tissue in real time. Finally, the real-time angular velocity and acceleration data are input into the nonlinear viscoelastic prediction model for iterative calculation to obtain stress-strain relationship data including instantaneous stress response and delayed strain components. This iterative calculation method based on dynamic data can update the prediction results of the model in real time to ensure its high consistency with the actual physiological state. The acquisition of instantaneous stress response and delayed strain component data provides an accurate mechanical basis for the subsequent lung tissue deformation analysis and puncture path planning, so that the puncture path can accurately adapt to the dynamic changes of lung tissue. Through the above steps, this method successfully constructs a nonlinear viscoelastic model that can accurately predict the dynamic deformation of lung tissue. This model not only fully considers the complex mechanical properties of lung tissue, but also combines real-time respiratory movement data to provide scientific and precise mechanical support for path planning of lung puncture surgery, significantly improving the safety and success rate of the puncture operation.

[0064] Preferably, step S25 comprises the following steps:

[0065] Step S251: inputting the real-time angular velocity data into the nonlinear viscoelastic prediction model of lung tissue, and calculating the instantaneous stress response data through the relationship between the angular velocity and the rotational deformation of the tissue;

[0066] In the embodiment of the present invention, in this step, the real-time angular velocity data is captured by infrared marker points and collected by a high-speed infrared imaging system at a sampling frequency of 100 Hz. The angular velocity data is substituted as an input parameter into the lung tissue non-linear viscoelasticity prediction model, and the model contains a coupling relationship formula between the angular velocity and the rotational deformation: Wherein, is the instantaneous stress response, is the moment of inertia of the tissue, is the angular velocity. The standard lung tissue moment of inertia is measured to be 0.0125 kg·m² through a calibration experiment, and the Python script is used to calculate the instantaneous stress response in real time. Taking a certain test data as an example, when the angular velocity fluctuates in the range of 0.5 rad / s to 1.5 rad / s, the calculated instantaneous stress response data shows a linear change trend, and the maximum value is about 0.01875 Pa.

[0067] Step S252: Input the real-time acceleration data into the lung tissue non-linear viscoelasticity prediction model, and calculate the delayed strain component data according to the relationship between the acceleration and the tissue inertial stress;

[0068] In the embodiment of the present invention, the real-time acceleration data is from the abdominal marker points and obtained by differentiating the displacement related to time. The data sampling rate is the same as that of the angular velocity, which is 100 Hz. Input the acceleration data into the lung tissue non-linear viscoelasticity prediction model, and use the following formula to calculate the inertial stress:

[0069] , wherein, is the inertial stress, is the tissue density, is the acceleration. The standard lung tissue density measured in the laboratory is 1.04 g / cm³ (about 1040 kg / m³). The acceleration data is processed by writing a MATLAB script to calculate the delayed strain component. Taking the actual acceleration data as an example, when the acceleration peak value is 0.3 m / s², the inertial stress is about 0.312 Pa and gradually decays to 0.05 Pa over time.

[0070] Step S253: Calibrate the instantaneous stress response data and the delayed strain component data to obtain the initial stress-strain relationship data;

[0071] In the embodiment of the present invention, the instantaneous stress response data in step S251 and the delayed strain component data in step S252 are merged, and the time alignment of the two sets of data is performed by an interpolation method. Linear interpolation is used to ensure that the time step is consistent (0.01 s). Subsequently, the merged data is calibrated using the standard lung tissue model data. For example, normalization processing is adopted to ensure that the data range is between [0, 1]. The calibrated initial stress-strain relationship data is verified by plotting the stress-strain curve. The curve shows linearity in the initial stage and then gradually transitions to non-linear characteristics, which conforms to the expected behavior of the non-linear viscoelastic model.

[0072] Step S254: Construct an error evaluation function based on the initial stress-strain relationship data to obtain model prediction error data;

[0073] In the embodiment of the present invention, an error evaluation function is constructed according to the initial stress-strain relationship data for calculating the difference between the model prediction value and the real experimental data. The form of the error evaluation function is: , where is the model predicted stress value, is the experimentally measured stress value, is the number of time steps. Batch data calculations are performed through MATLAB to obtain the model prediction error data. For example, the error value under a certain set of experimental conditions is 0.12 Pa, indicating that there is a certain room for optimization of the model.

[0074] Step S255: Compare the model prediction error data with a preset threshold. When the error is greater than the threshold, optimize and adjust the parameters of the lung tissue non-linear viscoelastic prediction model and return to step S23 for recalculation until the error is less than the threshold to obtain the final stress-strain relationship data.

[0075] In the embodiment of the present invention, the error data output by the error evaluation function is compared with a preset threshold (such as 0.05 Pa). When the error is greater than the threshold, the parameters in the lung tissue non-linear viscoelastic prediction model are adjusted, such as the elastic modulus and viscosity coefficient, and the calculation process of steps S23 to S254 is re-executed. The genetic algorithm (GA) is used in the optimization process, with the initial population size set to 50, the maximum number of iterations set to 200, and the goal is to minimize the error function. After 3 rounds of optimization iterations, the error is reduced to 0.03 Pa. The final obtained stress-strain relationship data is verified by experiments, and the error range is less than 2%. The optimized data curve is highly consistent with the real experimental curve and can be used for actual lung tissue dynamic analysis applications.

[0076] The present invention inputs real-time angular velocity data into a non-linear viscoelastic prediction model, and calculates instantaneous stress response data based on the relationship between angular velocity and tissue rotational deformation. This process can accurately capture the instantaneous mechanical response of lung tissue during respiratory movement, reflecting the instantaneous stress state of the tissue during dynamic deformation. By introducing angular velocity data, the model can more accurately simulate the rotational deformation characteristics of lung tissue, thereby providing key mechanical parameters for subsequent deformation prediction. Real-time acceleration data is input into the model, and by analyzing the relationship between acceleration and tissue inertial stress, the delayed strain component data is calculated. This process further improves the model's description of the dynamic behavior of lung tissue, especially considering the delayed deformation characteristics of the tissue under inertial action. The introduction of the delayed strain component enables the model to more comprehensively reflect the complex mechanical behavior of lung tissue during the respiratory cycle, providing more abundant information for deformation analysis. The instantaneous stress response data and the delayed strain component data are calibrated to obtain the initial stress-strain relationship data. This calibration process ensures the accuracy and consistency of the model output data, eliminating possible systematic errors or data biases. Through data calibration, the model can more realistically reflect the actual mechanical behavior of lung tissue, providing a reliable basis for subsequent error evaluation and optimization. An error evaluation function is constructed based on the initial stress-strain relationship data, and the model prediction error is calculated. This process provides a quantitative basis for the optimization of the model, enabling the accuracy of the model to be intuitively reflected by the error data. The introduction of the error evaluation function enables the model to self-evaluate the accuracy of its prediction results, providing a clear direction for subsequent parameter optimization. The model prediction error is compared with a preset threshold. When the error exceeds the threshold, the model parameters are optimized and adjusted, and the process returns to step S23 for recalculation until the error is less than the threshold. This iterative optimization mechanism ensures that the final output result of the model has high precision and high reliability. By dynamically adjusting the model parameters and repeatedly verifying, the model can continuously approximate the actual mechanical behavior of lung tissue, thereby providing accurate deformation prediction data for puncture path planning. Through the above steps, the method can not only accurately simulate the dynamic mechanical behavior of lung tissue during respiration, but also ensure the high precision and reliability of the model through a strict error evaluation and optimization mechanism. This refined model optimization process provides scientific and accurate mechanical support for the path planning of lung puncture surgery, significantly improving the safety and success rate of the surgery.

[0077] Preferably, step S3 includes the following steps:

[0078] Step S31: Obtain the coordinate data of the puncture target point and the size range data of the puncture target area;

[0079] In the embodiments of the present invention, the chest CT scan data of a patient is processed by a medical imaging workstation (such as 3D Slicer) to mark the puncture target points. The doctor manually selects the target points according to the anatomical features of the target area (such as the tumor center or the lesion area) and records their three-dimensional coordinate data, for example, (50.2mm, 32.8mm, 78.5mm). The data on the size range of the target area is automatically calculated by an image segmentation algorithm. Based on the segmented tumor volume, the maximum diameter and the minimum diameter of the target area are obtained. For example, the size range of the puncture target area of a certain patient is from 15mm to 22mm in diameter, and the data is stored in JSON format for subsequent processing.

[0080] Step S32: Determine a spherical region with a radius of 20mm centered on the target point according to the coordinate data of the puncture target point, so as to obtain the data on the range of the central region;

[0081] In the embodiments of the present invention, taking the target point coordinates obtained in step S31 as the center of the sphere and setting the sphere radius to 20mm, a three-dimensional spherical region is generated using the sphere geometric formula: , a spherical grid model is generated by calling the NumPy and Matplotlib libraries through a Python script, with the target point (50.2, 32.8, 78.5) as the center of the sphere and the radius set to 20mm. The data structure of the spherical region is stored in the form of a coordinate set of grid points, containing 10,000 evenly distributed points. For example, in a certain simulation, the data on the range of the central region includes the coordinate ranges (30.2, 70.2), (12.8, 52.8), (58.5, 98.5).

[0082] Step S33: Extract the chest wall boundary according to the CT scan sequence data, so as to obtain the data on the range of the peripheral region within 30mm from the inner surface of the chest wall;

[0083] In the embodiments of the present invention, after loading the CT scan data, a gradient-based boundary extraction algorithm (such as the Sobel operator) is used to extract the inner surface of the chest wall. The chest wall boundary points identified in each layer of the CT image are combined into a three-dimensional surface model, and the region within 30mm is expanded through a dilation operation. The dilation operation is implemented using the scikit-image library in Python, and the expanded range of the peripheral region is described by point cloud data. For example, in a certain extraction process, the peripheral region contains a total of 45,000 points and is stored in PLY format for subsequent stress analysis.

[0084] Step S34: Calculate the stress distribution through high-density and low-density grid division according to the data on the range of the central region and the data on the range of the peripheral region respectively, so as to obtain the first deformation field data of the central region and the second deformation field data of the peripheral region;

[0085] In the embodiment of the present invention, since the central region is close to the target point and requires high-precision calculation, a grid spacing of 2 mm is used for division; a low-density division is performed on the peripheral region with a grid spacing of 5 mm. Two sets of grid data are loaded through a finite element analysis software (such as ANSYS or COMSOL), and the elastic modulus (10 kPa) and Poisson's ratio (0.45) of the tissue are set as material parameters. When the simulated loading pressure is 100 Pa, the first deformation field data of the central region and the second deformation field data of the peripheral region are calculated respectively through a stress field solver. The results show that the maximum displacement of the central region is 0.8 mm, while the maximum displacement of the peripheral region is 0.3 mm. The two sets of deformation field data are exported in VTK format for subsequent processing.

[0086] Step S35: Refine the regions with large stress gradients according to the first deformation field data of the central region and the second deformation field data of the peripheral region through adaptive optimized meshing, and finally calculate the third deformation field data of the transition region.

[0087] In the embodiment of the present invention, the two sets of deformation field data in step S34 are input into the grid optimization module, and an adaptive algorithm is used to identify the regions where the stress gradient change exceeds 5 Pa / mm. In these regions, the grid spacing is further refined to 1 mm to enhance the calculation accuracy. Subsequently, the deformation field is recalculated through a finite element solver to generate the third deformation field data of the transition region. Through result analysis, it is found that the stress distribution in the transition region is highly continuous with that in the central region and the peripheral region, indicating that the optimized grid can accurately capture the stress change characteristics. The final third deformation field data of the transition region is visually displayed through ParaView, and the accuracy meets the requirements of medical applications and is used for subsequent puncture path planning and risk assessment.

[0088] The present invention obtains the coordinate data of the puncture target point and the size range data of the puncture target area, providing a clear positioning basis for subsequent regional division and deformation field calculation. This process ensures that the puncture path planning can accurately focus on the target area while considering the overall range of the target area, laying a foundation for subsequent mesh division and deformation analysis. A spherical central area with a radius of 20 mm is defined centered on the puncture target point, and this area division fully considers the accuracy requirements of the puncture operation. By defining the central area range, resources can be concentrated to perform high-precision modeling and deformation analysis on the key puncture area, ensuring the accuracy and safety of the puncture path. The chest wall boundary is extracted from the CT scan sequence data, and the peripheral area range within 30 mm from the inner surface of the chest wall is determined. This process not only considers the relative position of the puncture path and the chest wall but also provides boundary conditions for subsequent deformation field calculation. The division of the peripheral area can effectively simulate the mechanical environment around the puncture path and provide comprehensive mechanical support for the puncture path planning. High-density and low-density mesh divisions are respectively adopted for the central area and the peripheral area, and the stress distribution is calculated. This zoning mesh strategy can efficiently balance the calculation accuracy and resource consumption. The high-density mesh division of the central area ensures high-precision deformation analysis of the key puncture area, while the low-density mesh division of the peripheral area reduces the calculation burden and retains the characteristics of the overall mechanical environment. Through this strategy, the first deformation field data of the central area and the second deformation field data of the peripheral area can be obtained quickly and accurately. Finally, the transition area is refined by adaptively optimizing the mesh, further improving the calculation of the deformation field. Adaptive mesh division can dynamically adjust the mesh density according to the change of the stress gradient, ensuring refined modeling of areas with large stress gradients. This dynamic optimization strategy not only improves the accuracy of the deformation field calculation but also enhances the adaptability of the model to complex mechanical environments, and finally obtains the third deformation field data of the transition area. Through the above steps, the method realizes the refined deformation field calculation of the key puncture area, the peripheral area, and the transition area. This zoning strategy and adaptive optimization mechanism not only improve the efficiency and accuracy of the deformation analysis but also provide comprehensive and accurate mechanical support for the puncture path planning, significantly enhancing the accuracy and safety of the puncture surgery.

[0089] Preferably, step S34 includes the following steps:

[0090] Step S341: Perform a 0.5-mm high-density mesh division on the central area range data to obtain the initial central area mesh data;

[0091] In the embodiment of the present invention, for the central region range data generated in step S32, a spherical region with the target point as the center of the sphere and a radius of 20 mm is used as the computational domain, and high-density mesh division is carried out through a mesh generation tool (such as GMSH or ANSYS Meshing). The mesh spacing is set to 0.5 mm, and the tetrahedral element division method is adopted to ensure that the mesh can carefully capture the possible stress gradient changes within the region. The specific operations include importing the geometric data of the spherical region, setting the division parameters (element size is 0.5 mm), and enabling the automatic smoothing function to optimize the mesh quality. The generated initial central region mesh data contains approximately 65,000 mesh elements, and the shape quality index (Skewness) of each element is controlled below 0.2 to ensure the calculation accuracy. In a typical simulation, the generated mesh is stored in VTK format for subsequent calculations.

[0092] Step S342: Calculate the stress distribution in the central region based on the initial central region mesh data and the stress-strain relationship data, so as to obtain the first deformation field data of the central region;

[0093] In the embodiment of the present invention, the initial central region mesh data generated in step S341 is imported into a finite element analysis software (such as COMSOL Multiphysics). In the material property settings, according to the mechanical properties of the lung tissue, the elastic modulus is defined as 10 kPa, the Poisson's ratio is 0.45, and the nonlinear viscoelastic model is enabled. The stress-strain relationship data is loaded, and the boundary condition of applying an external uniform pressure of 150 Pa on the surface of the spherical region is set. Through the built-in solver of the software (such as Direct SparseSolver), a static analysis is carried out to calculate the stress distribution and deformation response within the region. The output first deformation field data of the central region includes the stress values and displacement vectors of each mesh point. For example, the maximum deformation at the target point is 0.82 mm, and the results are saved in CSV and VTK formats for subsequent analysis.

[0094] Step S343: Perform a low-density mesh division with a 2-mm interval on the peripheral region range data to obtain the initial peripheral region mesh data;

[0095] In the embodiment of the present invention, according to the peripheral region range data generated in step S33, a low-density grid division is performed on the region within 30 mm from the inner surface of the chest wall. Using the same grid division tool, the grid spacing is set to 2 mm, and the division strategy adopts a mixed tetrahedral element division to ensure the balance between calculation efficiency and result accuracy. The operation steps include: importing the geometric data of the peripheral region, setting the target element size to 2 mm, and enabling the boundary layer refinement function to optimize the grid quality at the region boundary. The generated initial peripheral region grid data contains approximately 12,000 grid elements, and the quality index (Aspect Ratio) of each element is less than 3, ensuring the rationality of the element shape. The finally generated grid is stored in the ABAQUS input file format (.inp) to facilitate compatibility with other calculation software.

[0096] Step S344: Calculate the stress distribution in the peripheral region according to the initial peripheral region grid data and the stress-strain relationship data, so as to obtain the second deformation field data of the peripheral region.

[0097] In the embodiment of the present invention, the peripheral region grid data generated in step S343 is imported into the finite element analysis software. Consistent with step S342, the lung tissue material parameters and the stress-strain relationship model are set. In terms of boundary conditions, a simulated internal physiological pressure gradient is applied to the peripheral region, and the pressure gradually increases from 80 Pa to 120 Pa from the chest wall boundary to the target region, simulating the real pressure distribution state within the simulation region. The model is calculated using a nonlinear static solver, and the stress and displacement data of each grid point are output. The second deformation field data of the peripheral region shows that the maximum displacement in this region is 0.35 mm, and the displacement and stress distributions are relatively uniform, meeting the requirements of subsequent analysis. The calculation results are exported in the HDF5 format and visually inspected through ParaView to verify the rationality and consistency of the distribution.

[0098] The present invention performs a high-density grid division with a 0.5-mm grid interval on the data of the central region. This high-precision grid division strategy can fully capture the subtle deformation characteristics near the puncture target point. Since the precision requirements for puncture operations are extremely high, the high-density grid division ensures accurate modeling of the mechanical behavior of the key region and provides high-quality basic data for subsequent deformation field calculations. This high-density grid division not only improves the accuracy of deformation analysis but also effectively reflects the complex mechanical behavior of the tissue around the puncture path. Based on the initial grid data of the central region and the stress-strain relationship data, the stress distribution is calculated to obtain the first deformation field data of the central region. This process, by combining a high-precision grid and an accurate mechanical model, can describe in detail the stress state and deformation characteristics around the puncture target point. This high-precision deformation field data provides a key basis for optimizing the puncture path and ensures that the puncture operation can be carried out on the safest and most effective path. A low-density grid division with a 2-mm grid interval is performed on the data of the peripheral region. This strategy takes into account the relatively lower importance of the peripheral region, but still needs to retain the characteristics of the overall mechanical environment. The low-density grid division can efficiently simulate the mechanical behavior of the peripheral region while reducing the consumption of computing resources, providing a macroscopic mechanical background for puncture path planning. This zoned grid strategy effectively balances the calculation accuracy and resource utilization efficiency. Based on the initial grid data of the peripheral region and the stress-strain relationship data, the stress distribution of the peripheral region is calculated to obtain the second deformation field data of the peripheral region. This process further improves the calculation of the overall deformation field, enabling the puncture path planning to fully consider the influence of the peripheral region on the puncture operation. By combining the deformation field data of the central region and the peripheral region, the safety and feasibility of the puncture path can be evaluated more comprehensively. Through the above steps, the method realizes the refined deformation field calculation of the key puncture region and the peripheral region. The high-density grid division ensures accurate modeling near the puncture target point, while the low-density grid division improves the calculation efficiency while retaining the overall mechanical characteristics.

[0099] Preferably, step S35 includes the following steps:

[0100] Step S351: Obtain the data of the transition region range other than the central region and the peripheral region according to the first deformation field data of the central region and the second deformation field data of the peripheral region;

[0101] In the embodiment of the present invention, the first deformation field data of the central region generated in step S342 and the second deformation field data of the peripheral region generated in step S344 are used to extract the range data of the transition region through the change of the deformation field gradient. The specific operation is to import the deformation field data of the central region and the peripheral region into the MATLAB environment, calculate the spatial gradient of the deformation field data using a gradient function (such as gradient), screen out the regions where the gradient value is between 0.05 and 0.2, and mark them as the transition regions. Subsequently, through Boolean operations, the parts belonging to the central region and the peripheral region are removed, and only the range data of these transition regions is retained. The output range data of the transition region is saved in the point cloud format and stored in the PLY format for subsequent mesh generation.

[0102] Step S352: Perform an initial mesh generation with a density of 1 mm on the range data of the transition region to obtain the initial mesh data of the transition region;

[0103] In the embodiment of the present invention, the range data of the transition region generated in step S351 is imported into a mesh generation software (such as GMSH or ABAQUS), and mesh generation is performed with an element spacing of 1 mm. The specific operations include: importing the point cloud data of the transition region, generating tetrahedral meshes through the Delaunay triangulation method, setting the target element size to 1 mm, and enabling the boundary smoothing function to optimize the mesh boundary quality. The generated initial mesh data contains approximately 25,000 elements, and the shape quality of each element is controlled below 0.3 to ensure the uniformity and stability of the mesh. The initial mesh data of the transition region is stored in the.inp format for subsequent stress distribution calculation.

[0104] Step S353: Calculate the local stress gradient data based on the initial mesh data of the transition region and the stress-strain relationship data;

[0105] In the embodiment of the present invention, the initial mesh data of the transition region generated in step S352 is imported into a finite element analysis software (such as ANSYS or COMSOL), and the stress-strain relationship data used in steps S342 and S344 is loaded. By setting the linear pressure distribution at the boundary of the transition region (from 120 Pa to 150 Pa), the stress response characteristics at different positions are simulated. A static solver is used to calculate the stress values and displacement values of each mesh point, and a gradient calculation tool (such as the Compute Gradient module) is used to extract the local stress gradient data. The calculation results show that the maximum value of the stress gradient in the transition region is approximately 0.18, and the minimum value is approximately 0.03. The local stress gradient data is exported in the CSV format for mesh optimization.

[0106] Step S354: Based on the local stress gradient data, adaptively optimize the transition region meshes in the initial transition region mesh data. When the stress gradient is greater than the preset threshold, subdivide the meshes to 0.5 mm, thereby obtaining the optimized mesh data for the transition region;

[0107] In the embodiment of the present invention, the local stress gradient data generated in step S353 is compared with a preset threshold (0.1), and the mesh elements with stress gradients greater than the threshold are screened out. For these high-gradient regions, an adaptive mesh optimization method is adopted to refine the element spacing from 1 mm to 0.5 mm, and the local remeshing technology (such as Adaptive Mesh Refinement) is used to achieve mesh refinement. During the meshing process, the boundary smoothness and element quality are ensured first, and mesh distortion caused by excessive refinement is avoided. The finally generated optimized mesh data contains approximately 40,000 elements, the element size in the high-gradient region is 0.5 mm, and the optimized mesh data is stored in the HDF5 format for subsequent stress calculation.

[0108] Step S355: Calculate the stress distribution in the transition region according to the optimized mesh data and the stress-strain relationship data, thereby obtaining the third deformation field data of the transition region.

[0109] In the embodiment of the present invention, the optimized mesh data generated in step S354 is imported into the finite element analysis software, the stress-strain relationship model and boundary conditions (the pressure gradually transitions from 120 Pa to 150 Pa) are loaded, and the stress distribution in the transition region is calculated. A nonlinear static solver is used to accurately solve the stress and deformation at each mesh point, especially paying attention to the stress concentration phenomenon in the high-gradient region. The calculation results show that the maximum displacement in the transition region is 0.6 mm, the minimum displacement is 0.1 mm, and the deformation field distribution is uniform and continuous. The output third deformation field data of the transition region is stored in the VTK format and verified visually through ParaView to ensure the accuracy and consistency of the results.

[0110] The present invention determines the range of the transition region by analyzing the deformation field data of the central region and the peripheral region. This process ensures that the definition of the transition region matches the deformation characteristics of the central and peripheral regions, providing clear boundary conditions for subsequent mesh generation and stress distribution calculation. The precise division of the transition region makes the modeling of the deformation field more complete and can better reflect the mechanical transition characteristics of the lung tissue between different regions. An initial mesh with a density of 1 mm is generated for the transition region, providing a basic framework for the deformation analysis of the transition region. This moderate mesh density can capture the mechanical behavior of the transition region initially while ensuring computational efficiency. The establishment of the initial mesh provides the necessary data structure for subsequent stress gradient calculation and mesh optimization. By combining the initial transition region mesh data and stress-strain relationship data, the local stress gradient is calculated. This process can identify the key regions with significant stress changes in the transition region, providing a quantitative basis for mesh optimization. The calculation of the local stress gradient enables the model to respond dynamically to the mechanical complexity of the transition region, ensuring that subsequent mesh optimization can accurately focus on the high stress gradient regions. The transition region mesh is adaptively optimized based on the local stress gradient data. When the stress gradient is greater than the preset threshold, the mesh is further refined to 0.5 mm, resulting in optimized mesh data. This adaptive optimization strategy can dynamically adjust the mesh density, ensuring high-precision modeling in the high stress gradient regions of the transition region while avoiding excessive consumption of computational resources in the low stress gradient regions. Through this optimization, the deformation field modeling of the transition region can better adapt to the complex mechanical environment. According to the optimized mesh data and stress-strain relationship data, the stress distribution in the transition region is calculated to obtain the third deformation field data. This process further improves the deformation field modeling of the transition region, making the deformation field of the entire lung tissue more continuous and accurate. The generation of the third deformation field data for the transition region enables the puncture path planning to fully consider the mechanical characteristics of the transition region, thereby achieving a safer and more accurate puncture path design. Through the above steps, the method realizes the refined modeling of the deformation field of the transition region. From the precise division of the transition region to the establishment of the initial mesh, then to the adaptive mesh optimization based on the stress gradient, and finally to the stress distribution calculation, this series of processes not only improves the modeling accuracy of the deformation field of the transition region but also enhances the integrity and continuity of the deformation field of the entire lung tissue.

[0111] Preferably, step S4 includes the following steps:

[0112] Step S41: Calculate the gray gradient of adjacent images in the CT scan sequence data to obtain image gray gradient data, where the image gray gradient data includes x-direction gradient data, y-direction gradient data, and time-direction gradient data;

[0113] In the embodiment of the present invention, in this step, CT scan sequence data is imported (for example, 50 consecutive frames of chest CT images with a resolution of 512×512 pixels and an interval time of 0.5 seconds). The Sobel operator is used to calculate the gradients in the x - direction and y - direction of each frame of the image respectively to extract the spatial gradient information. At the same time, by using the change of the image sequence in the time dimension, the gray - level gradient in the time direction is calculated through the difference method (such as the backward difference formula). The calculation formulas are as follows: Gradient in the x - direction: , and convolution calculation is performed using a 3×3 Sobel operator; Gradient in the direction: , and convolution calculation is performed using the same method; Gradient in the time direction: , which is calculated through the backward difference formula, that is . The gradient calculation results are saved in the 3D matrix format (with a size of 512×512×50) for constructing the brightness constraint condition in the subsequent optical flow method.

[0114] Step S42: Construct initial motion vector field data based on the first deformation field data in the central region, the second deformation field data in the peripheral region, and the third deformation field data in the transition region;

[0115] In the embodiment of the present invention, the deformation field data of the central region, the peripheral region, and the transition region generated in steps S342, S344, and S355 are loaded into the MATLAB environment. Three - dimensional interpolation processing is performed on the deformation fields of different regions to unify the resolution to 1mm³ to ensure the spatial consistency of the data. The data of different regions are combined through a linear weighting method: ; where , and are the weight coefficients of the central region, the peripheral region, and the transition region respectively (the initial values are 0.4, 0.3, 0.3, which are set according to the importance of the region), is the deformation field data of each region. The calculation result is the initial motion vector field data, which is stored in the vector field format (such as a.vtk file).

[0116] Step S43: Establish a brightness - constancy constraint equation based on the image gray - level gradient data to obtain the brightness constraint condition data;

[0117] In the embodiment of the present invention, the gray - level gradient data calculated in step S41 is substituted into the brightness - constancy hypothesis model: , where , and are the components of the gray - level gradient in the direction, direction, and the time direction respectively, , is the motion vector field in the and Components in the

[0118] Step S44: Establish a velocity smoothing constraint equation based on the initial motion vector field data, thereby obtaining velocity constraint condition data;

[0119] In the embodiment of the present invention, the initial motion vector field data generated in step S42 is imported, and the velocity gradient is calculated for each grid cell: In a finite element analysis tool (such as COMSOL), the above equation is discretized into a system of linear equations, and an iterative solution method is used to generate velocity constraint condition data, and the result is saved in the.mat format.

[0120] Step S45: Combine the luminance constraint condition data and the velocity constraint condition data with weights to obtain an optical flow energy functional expression; and perform Euler-Lagrange equation derivation on the optical flow energy functional expression to obtain a control system of equations for the variational problem;

[0121] In the embodiment of the present invention, the energy functional of the optical flow method is constructed by combining the constraint conditions of steps S43 and S44:

[0122] ; where and are the weights of the luminance constraint and the velocity smoothing constraint (the initial values are 0.7 and 0.3), respectively. Through variational method, Euler-Lagrange derivation is performed on the energy functional to obtain the following control system of equations: Use a symbolic calculation tool (such as the syms module of MATLAB) to expand and simplify the equations, generate the final control system of equations, and output it in Latex format for verification.

[0123] Step S46: Numerically solve the control system of equations to obtain the optimized motion vector field data for each region;

[0124] In the embodiment of the present invention, the control system of equations derived in step S45 is imported into a solution tool (such as ANSYS or a custom Python script), boundary conditions and initial values are set, and the finite difference method (FDM) or the finite element method (FEM) is used for numerical solution. During the calculation process, the Gauss-Seidel iteration method is used to improve the convergence speed, and the initial iteration accuracy is set to . The solution result is the optimized motion vector field data, which is saved in the.hdf5 format.

[0125] Step S47: Perform data fusion based on the optimized motion vector field data of each region, eliminate the discontinuity between regions through boundary smoothing processing, construct an error function, and compare it with a preset threshold. When the error is greater than the threshold, return to step S45 to adjust the weight coefficient and recalculate until the error is less than the threshold, so as to obtain the tissue deformation field data.

[0126] In the embodiment of the present invention, the optimized motion vector field data of each region is imported into ParaView for data visualization to check the continuity of the boundary region. Through the boundary smoothing algorithm based on bicubic interpolation, the mutation of the vector field at the boundaries of different regions is eliminated. Then, define the error function: Compare the calculated error value with a preset threshold (for example, 0.001). If the error is greater than the threshold, return to step S45 to adjust and the ratio of, and recalculate until the accuracy requirement is met. The finally generated tissue deformation field data is saved in VTK format, and the continuity and accuracy of the three-dimensional deformation field are verified through ParaView.

[0127] The present invention obtains gradient data in the x-direction, y-direction, and time-direction by calculating the gray-scale gradients of adjacent images in a CT scan sequence. This process provides basic information for subsequent calculation of the motion vector field and can effectively capture the dynamic changes of lung tissue during respiratory motion. The introduction of the gray-scale gradient data enables the calculation of the motion vector field to be based on the brightness changes of the images, thereby more accurately reflecting the motion characteristics of the tissue. An initial motion vector field is constructed based on the deformation field data of the central region, peripheral region, and transitional region. This process integrates the deformation information of different regions and provides comprehensive initial conditions for subsequent optimization calculations. By combining the deformation field data of multiple regions, the overall motion state of the lung tissue can be more completely described, providing richer dynamic information for puncture path planning. A brightness constancy constraint equation and a velocity smoothness constraint equation are established respectively. The brightness constancy constraint uses the gray-scale gradient data to ensure that the calculation of the motion vector field conforms to the law of image brightness changes, while the velocity smoothness constraint guarantees the smoothness and continuity of the motion vector field. The introduction of these two constraint conditions makes the calculation of the motion vector field more stable and reliable, and can effectively avoid calculation errors caused by noise or local anomalies. By weighted combination of the brightness constraint condition and the velocity constraint condition, an optical flow method energy functional expression is constructed, and the control equations of the variational problem are derived. This process combines the deformation field data of multiple regions with the image brightness changes to form a unified optimization framework. Through the derivation of the Euler-Lagrange equation, it can be theoretically guaranteed that the calculation result of the motion vector field has an optimal solution, thereby improving the accuracy of the deformation field modeling. By numerically solving the control equations, the optimized motion vector field data of each region are obtained. This process further refines the calculation of the motion vector field and ensures the accuracy and consistency of the motion vectors of each region. The optimized motion vector field can more realistically reflect the dynamic changes of the lung tissue during respiration and provides reliable dynamic data support for puncture path planning. Through data fusion and boundary smoothing processing, the discontinuity between different regions is eliminated, and the motion vector field is iteratively optimized by comparing the error function with a preset threshold. This process not only improves the overall quality of the deformation field data but also ensures the accuracy of the final result by dynamically adjusting the weight coefficients. The optimized tissue deformation field data can more accurately reflect the dynamic behavior of the lung tissue and provide high-quality input for puncture path planning. Through the above steps, the method realizes the accurate calculation of the motion vector field of the lung tissue and the optimized modeling of the deformation field. From the acquisition of the gray-scale gradient data to the integration of the deformation field data of multiple regions, and then to the establishment of the constraint conditions and the optimized calculation, this series of processes not only improves the accuracy and stability of the deformation field modeling but also provides comprehensive and reliable dynamic data support for puncture path planning, significantly improving the accuracy and safety of the puncture surgery.

[0128] Preferably, step S5 includes the following steps:

[0129] Step S51: Perform gray value threshold segmentation on the CT sequence data to obtain initial lung contour data; extract lung surface feature line data according to the initial lung contour data, where the lung surface feature line data includes the surface projections of the pulmonary fissures and the main blood vessel courses.

[0130] In the embodiment of the present invention, gray value threshold segmentation is performed on the CT sequence data, and the Otsu adaptive threshold segmentation method is used to separate the lung region from other tissues. Specifically, the optimal segmentation threshold T is calculated for the gray histogram in each CT image, and the lung region (low-density region) is extracted through T to obtain the initial lung contour data. Subsequently, the contour tracking algorithm is used to extract the outer boundary of the lung, and the lung surface feature line data is extracted by combining the curvature-based feature extraction method. In actual operation, by comparing the gradient changes in the pulmonary fissure and main blood vessel regions, the Canny edge detection algorithm is used to further identify the surface projections of the pulmonary fissures and the main blood vessels, and high-precision lung surface feature line data is generated. During the processing, the resolution of the CT image is 512×512, and the pixel pitch is 0.5 mm.

[0131] Step S52: Perform time series analysis on the tissue deformation field data to obtain tissue deformation law data at different respiratory phases; establish an elastic deformation model based on the water bubble animation simulation algorithm, regarding the tissue as a droplet with surface tension, so as to obtain surface tension parameter data.

[0132] In the embodiment of the present invention, the tissue deformation field data obtained in step S47 is input into the time series analysis model, and the Fourier transform is used to extract the frequency characteristics and amplitude characteristics of the tissue at different respiratory phases to obtain the tissue deformation law data. On this basis, combined with the water bubble animation simulation algorithm, assuming that the tissue has elastic deformation characteristics similar to a droplet, by establishing a relational formula for the surface tension and volume change , the surface tension parameter data at different respiratory phases is calculated. Among them, is the surface tension coefficient, is the surface area change amount, is the volume change amount. For lung tissue, the value range of the surface tension coefficient is 0.02–0.04 N / m.

[0133] Step S53: Construct a dynamic equation describing the deformation of lung tissue according to the surface tension parameter data and the tissue deformation law data, so as to obtain deformation dynamics model data.

[0134] In the embodiment of the present invention, according to the surface tension parameter data and the tissue deformation law data obtained in step S52, the finite element method is used to construct the dynamic equation of lung tissue deformation. The dynamic equation is based on the stress balance formula where is the stress tensor, is the external force, is the density, and is the acceleration. Combining the elastic modulus of lung tissue (about 5–10 kPa) and the Poisson's ratio (taking the value of 0.35), the dynamic equation is discretized to obtain the deformation dynamics model data.

[0135] Step S54: Numerically solve the deformation dynamics model data to obtain the three-dimensional grid node position data at different breathing phases; perform surface reconstruction based on the three-dimensional grid node position data and maintain the original grid density in the characteristic line region to obtain the initial three-dimensional model data;

[0136] In the embodiment of the present invention, the deformation dynamics model data is input into a numerical solver (such as COMSOL Multiphysics) to simulate the tissue deformation at different breathing phases, and the three-dimensional grid node position data at each moment is output. During the grid reconstruction process, the grid density in the characteristic line region of the lung surface is kept unchanged, and the initial three-dimensional model data is generated by the Marching Cubes algorithm. To ensure the accuracy of the model, the grid resolution in the characteristic line region is set to 0.2 mm, while that in other regions is 0.5 mm.

[0137] Step S55: Simplify the grid of the initial three-dimensional model data in the non-characteristic line region, reduce the number of grids to 50% of the original by the quadratic error metric criterion to obtain the optimized three-dimensional model data; identify the spatial distribution of the ribs and large blood vessels in the optimized three-dimensional model data to obtain the data of the prohibited penetration region;

[0138] In the embodiment of the present invention, the quadratic error metric criterion is used to simplify the grid of the initial three-dimensional model data in the non-characteristic line region, reduce the number of grids to 50% of the original, and retain the geometric accuracy of the key structures during the simplification process. The spatial distribution of the ribs and large blood vessels is identified by the principal curvature analysis method and used as the prohibited penetration region to generate the data of the prohibited penetration region. In actual processing, the boundary error of the prohibited penetration region is less than 0.5 mm to ensure the reliability of the safety assessment of the puncture path.

[0139] Step S56: Calculate the feasible puncture path based on the data of the prohibited penetration region and the breathing phase mapping relationship to obtain the data of the initial path set; perform a safety assessment on the data of the initial path set to obtain the path risk assessment data, where the safety assessment is specifically to calculate the minimum distance between each path and the prohibited penetration region at different breathing phases;

[0140] In an embodiment of the present invention, by combining the data of the forbidden penetration area and the deformation mapping relationship between different breathing phases, a shortest path calculation method based on the Dijkstra algorithm is adopted to analyze the spatial distance and deformation trajectory of each path, and initial path set data is generated. For each path, the minimum distance from the forbidden penetration area at each breathing phase is calculated, and the minimum distance is used as a safety evaluation index to generate path risk assessment data. The minimum safety distance threshold is set to 3 mm, and the evaluation results are used to screen high-risk paths.

[0141] Step S57: Screen and sort the initial path set according to the path risk assessment data, so as to obtain the final puncture path planning model data.

[0142] In an embodiment of the present invention, the initial path set is screened according to the path risk assessment data, and the paths with a safety distance less than 3 mm are excluded. Subsequently, the remaining paths are sorted according to the risk assessment scores, and the optimal path is selected as the final puncture path planning model data. During this process, the path connection points are smoothed by the boundary smooth interpolation method to ensure the continuity of the path and the accuracy of the puncture. The finally generated path planning data is used to guide the actual surgical operation and is clinically verified. The verification results show that the deviation of the path planning is within 0.5 mm.

[0143] The present invention processes CT sequence data through grayscale value threshold segmentation to extract initial lung contour and surface feature line data, including the surface projections of pulmonary fissures and major blood vessels. This process provides an accurate anatomical structure basis for subsequent 3D modeling, ensuring that the model can truly reflect the geometric features and important anatomical landmarks of the lungs. Perform temporal analysis on tissue deformation field data, and establish an elastic deformation model in combination with the blister animation simulation algorithm. This modeling method based on dynamic deformation laws can effectively capture the deformation characteristics of lung tissue at different respiratory phases. At the same time, by introducing surface tension parameters, it further improves the model's ability to describe tissue dynamic behavior. By constructing and solving the deformation dynamics model, the 3D mesh node position data at different respiratory phases is obtained, and surface reconstruction is performed. This process not only ensures the dynamics and accuracy of the 3D model, but also maintains the original mesh density in the feature line area, enabling the details of important anatomical structures to be retained. This refined modeling strategy provides a high-quality 3D anatomical background for puncture path planning. Simplify the 3D model of the non-feature line area to reduce the computational burden, and at the same time, identify the prohibited puncture areas through the spatial distribution of ribs and major blood vessels. This optimization process not only improves the computational efficiency of the model, but also provides key information for the safety assessment of the puncture path. Further, based on the data of the prohibited puncture area and the respiratory phase mapping relationship, calculate the feasible puncture paths, and perform safety assessment and screening on the path set. By calculating the minimum distance between each path and the prohibited puncture area, the risk of the path can be effectively evaluated, and finally a safe and feasible puncture path planning model is selected.

[0144] Preferably, the present invention further provides a 3D imaging reconstruction and positioning system for lung puncture, which is used to execute the above-mentioned 3D imaging reconstruction and positioning method for lung puncture. The 3D imaging reconstruction and positioning system for lung puncture includes:

[0145] A data acquisition module, configured to obtain the CT scan sequence data of the patient and the multi-point infrared marker data of the chest wall, wherein the multi-point infrared marker data of the chest wall includes the respiratory reference position data, the real-time angular velocity data of the chest wall movement state, and the real-time acceleration data of each point on the chest wall;

[0146] A tissue mechanics modeling module, configured to predict the non-linear viscoelasticity of lung tissue based on the Maxwell model according to the multi-point infrared marker data of the chest wall, and input the real-time angular velocity data and the real-time acceleration data into the non-linear viscoelasticity prediction model of lung tissue for iterative calculation to obtain the stress-strain relationship data of lung tissue, wherein the stress-strain relationship data includes the instantaneous stress response data and the delayed strain component data;

[0147] The partitioned grid division module is used to perform partitioned grid division on lung tissue based on stress-strain relationship data, obtaining first deformation field data of the central region based on high-density grids, second deformation field data of the peripheral region based on low-density grids, and third deformation field data of the transition region constructed with adaptive grids based on local stress gradients in the remaining region;

[0148] The motion vector calculation module is used to calculate the motion vector field for the first deformation field data of the central region, the second deformation field data of the peripheral region, and the third deformation field data of the transition region by using the optical flow method, construct an energy functional according to the brightness constancy constraint and the velocity smoothness constraint, and solve it by variational method to obtain the optimized tissue deformation field data;

[0149] The three-dimensional reconstruction and path planning module is used to apply the tissue deformation field data to CT scan sequence data and perform three-dimensional reconstruction based on the bubble animation simulation algorithm, simplify the grids in the non-feature line region while maintaining the grid density in the feature line region, and obtain the puncture path planning model data with respiratory phase mapping.

[0150] In the present invention, the data acquisition module provides high-quality basic information for subsequent modeling and analysis by obtaining the CT scan sequence data of the patient and the multi-point infrared marker data of the chest wall. The CT scan sequence data can record the anatomical structure of the lungs in detail, while the multi-point infrared marker data of the chest wall accurately captures the impact of respiratory movement on lung tissue through respiratory reference positions, real-time angular velocity, and acceleration information. The fusion of such multi-dimensional data enables the system to comprehensively consider the dynamic changes of lung tissue during the respiratory cycle, providing accurate real-time data support for puncture path planning. The tissue mechanics modeling module constructs a non-linear viscoelastic prediction model of lung tissue based on the Maxwell model and performs iterative calculations through real-time angular velocity and acceleration data. This process not only considers the complex mechanical properties of lung tissue but also optimizes the accuracy of the model through real-time feedback of dynamic data. The generation of instantaneous stress response and delayed strain component data provides key mechanical parameters for subsequent deformation field calculations, enabling the system to accurately simulate the dynamic deformation behavior of lung tissue during respiration and providing a reliable mechanical basis for puncture path planning. The partitioned mesh generation module performs partitioned mesh generation on lung tissue according to stress-strain relationship data, generating deformation field data for the central region, peripheral region, and transition region respectively. This partitioning strategy achieves a balance between computational accuracy and efficiency by focusing on key puncture regions with high-density meshes, covering peripheral regions with low-density meshes, and optimizing the transition region with adaptive meshes based on local stress gradients. Partitioned mesh generation not only improves the accuracy of deformation field calculations but also provides an efficient mesh basis for subsequent 3D reconstruction and path planning, ensuring the scientificity and feasibility of puncture path planning. The motion vector calculation module uses the optical flow method combined with brightness constancy constraints and velocity smoothness constraints to calculate the motion vector field for the deformation field data of each region. By constructing an energy functional and solving it using the variational method, the system can optimize the motion vector field and eliminate noise and discontinuities. The optimized tissue deformation field data can accurately reflect the dynamic changes of lung tissue during respiration, providing high-quality dynamic information support for puncture path planning and further enhancing the accuracy and safety of the puncture path. Finally, the 3D reconstruction and path planning module applies the tissue deformation field data to the CT scan sequence data and performs 3D reconstruction through the bubble animation simulation algorithm. During the reconstruction process, the system maintains the mesh density in the feature line region while simplifying the mesh in the non-feature line region, generating a puncture path planning model with respiratory phase mapping. This 3D reconstruction strategy not only retains the detailed information of important anatomical structures but also ensures that the puncture path can adapt to the dynamic changes of lung tissue through respiratory phase mapping. The finally generated puncture path planning model can provide intuitive and reliable navigation support for clinicians, significantly improving the success rate and safety of puncture surgery. Through the collaborative work of the above modules, the lung puncture 3D image reconstruction and positioning system realizes the full-process optimization from data acquisition to puncture path planning.Each module not only complements each other technically, but also forms an organic whole functionally, providing precise and efficient technical support for lung puncture surgery and promoting the development of clinical puncture technology to a higher level.

[0151] Therefore, from any perspective, the embodiments should be regarded as exemplary and non-limiting. The scope of the present invention is defined by the appended claims rather than the above description. Thus, all changes falling within the meaning and scope of the equivalent elements of the application document are intended to be encompassed within the present invention.

[0152] The above are only specific embodiments of the present invention, enabling those skilled in the art to understand or implement the present invention. Various modifications to these embodiments will be obvious to those skilled in the art. The general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to these embodiments shown herein, but rather to the broadest scope consistent with the principles and novel features invented herein.

Claims

1. A three-dimensional image reconstruction and localization method for lung puncture, characterized in that, It includes the following steps: Step S1: Obtain the CT scan sequence data of the patient and the multi-point infrared marker data of the chest wall, where the multi-point infrared marker data of the chest wall includes the respiratory reference position data, the real-time angular velocity data of the chest wall movement state, and the real-time acceleration data of each point on the chest wall; Step S2: Based on the multi-point infrared marker data of the chest wall, a non-linear viscoelastic prediction model of lung tissue based on the Maxwell model is used, and the real-time angular velocity data and the real-time acceleration data are input into the non-linear viscoelastic prediction model of lung tissue for iterative calculation to obtain the stress-strain relationship data of lung tissue, where the stress-strain relationship data includes the instantaneous stress response data and the delayed strain component data; Step S3: Based on the stress-strain relationship data, zoning meshing of the lung tissue is performed to obtain the first deformation field data of the central region based on high-density meshes, the second deformation field data of the peripheral region based on low-density meshes, and the third deformation field data of the transition region constructed by an adaptive mesh based on the local stress gradient in the remaining region; Step S4: Use the optical flow method to calculate the motion vector field of the first deformation field data of the central region, the second deformation field data of the peripheral region, and the third deformation field data of the transition region, and construct an energy functional according to the brightness constancy constraint and the velocity smoothness constraint, and solve it by the variational method to obtain the optimized tissue deformation field data; Step S5: Apply the tissue deformation field data to the CT scan sequence data and perform 3D reconstruction based on the bubble animation simulation algorithm, simplify the mesh in the non-feature line region while maintaining the mesh density in the feature line region, and obtain the puncture path planning model data with respiratory phase mapping.

2. The three-dimensional image reconstruction and positioning method for lung puncture according to claim 1, wherein, Step S1 includes the following steps: Step S11: Collect the CT scan sequence data during the patient's respiratory cycle; Step S12: Arrange multiple infrared reflection marker points on the surface of the patient's chest wall to obtain the infrared reflection marker point data, where a marker point is arranged at the manubrium sterni position to obtain the respiratory reference position information, marker points are arranged at the ends of the left and right seventh ribs respectively to obtain the angular velocity information, and a marker point is arranged 5 cm above the umbilicus of the abdomen to obtain the acceleration information; Step S13: Use a high-speed infrared imaging system to dynamically track the infrared reflection marker point data, and obtain the three-dimensional spatial coordinate data of each marker point through an image registration algorithm; Step S14: Smooth the motion trajectories of the marker points based on the three-dimensional spatial coordinate data to obtain the respiratory reference position data with noise effects eliminated; Step S15: Perform derivative operations on the motion trajectories of the marker points based on the three-dimensional spatial coordinate data to obtain the real-time angular velocity data and the real-time acceleration data respectively; Step S16: Align the respiratory reference position data, the real-time angular velocity data, and the real-time acceleration data in time series to obtain the multi-point infrared marker data of the chest wall.

3. The method for three-dimensional image reconstruction and positioning of lung puncture according to claim 2, wherein, Step S2 includes the following steps: Step S21: Obtain the elastic modulus parameter data and viscosity coefficient parameter data of the standard human lung tissue; Step S22: Establish an elastic element model in the Maxwell model based on the elastic modulus parameter data of the lung tissue to obtain initial elastic response data; establish a viscous element model in the Maxwell model based on the viscous coefficient parameter data of the lung tissue to obtain initial viscous response data; Step S23: Combine the initial elastic response data and the initial viscous response data in series to obtain basic Maxwell unit model data; Step S24: Establish a multi-dimensional state equation based on the basic Maxwell unit model data and the chest wall multi-point infrared marker data, and perform numerical discretization to obtain a lung tissue non-linear viscoelastic prediction model; Step S25: Input the real-time angular velocity data and the real-time acceleration data into the lung tissue non-linear viscoelastic prediction model for iterative calculation to obtain lung tissue stress-strain relationship data, where the stress-strain relationship data includes instantaneous stress response data and delayed strain component data.

4. The method for three-dimensional image reconstruction and localization of lung puncture according to claim 3, wherein, Step S25 includes the following steps: Step S251: Input the real-time angular velocity data into the lung tissue non-linear viscoelastic prediction model, and calculate the instantaneous stress response data through the relationship between the angular velocity and the tissue rotational deformation; Step S252: Input the real-time acceleration data into the lung tissue non-linear viscoelastic prediction model, and calculate the delayed strain component data according to the relationship between the acceleration and the tissue inertial stress; Step S253: Calibrate the instantaneous stress response data and the delayed strain component data to obtain initial stress-strain relationship data; Step S254: Construct an error evaluation function based on the initial stress-strain relationship data to obtain model prediction error data; Step S255: Compare the model prediction error data with a preset threshold. When the error is greater than the threshold, optimize and adjust the parameters of the lung tissue non-linear viscoelastic prediction model and return to Step S23 for recalculation until the error is less than the threshold to obtain the final stress-strain relationship data.

5. The method for three-dimensional image reconstruction and positioning of lung puncture according to claim 4, wherein Step S3 includes the following steps: Step S31: Obtain the puncture target point coordinate data and the puncture target area size range data; Step S32: Determine a spherical area with a radius of 20 mm centered on the target point according to the puncture target point coordinate data to obtain the central area range data; Step S33: Extract the chest wall boundary according to the CT scan sequence data to obtain the peripheral area range data within 30 mm from the inner surface of the chest wall; Step S34: Calculate the stress distribution through high-density and low-density grid division according to the central area range data and the peripheral area range data respectively to obtain the first deformation field data of the central area and the second deformation field data of the peripheral area; Step S35: Refine the area with a larger stress gradient according to the first deformation field data of the central area and the second deformation field data of the peripheral area through adaptive optimized grid, and finally calculate the third deformation field data of the transition area.

6. The method for three-dimensional image reconstruction and positioning of lung puncture according to claim 5, wherein, Step S34 includes the following steps: Step S341: Perform 0.5 mm high-density grid division on the central area range data to obtain initial central area grid data; Step S342: Calculate the stress distribution in the central region based on the initial central region grid data and the stress-strain relationship data, so as to obtain the first deformation field data of the central region; Step S343: Perform a 2-mm low-density grid division on the peripheral region range data to obtain the initial peripheral region grid data; Step S344: Calculate the stress distribution in the peripheral region based on the initial peripheral region grid data and the stress-strain relationship data, so as to obtain the second deformation field data of the peripheral region.

7. The method for three-dimensional image reconstruction and positioning of lung puncture according to claim 6, wherein Step S35 includes the following steps: Step S351: Obtain the transition region range data other than the central region and the peripheral region based on the first deformation field data of the central region and the second deformation field data of the peripheral region; Step S352: Perform an initial grid division with a density of 1 mm on the transition region range data to obtain the initial transition region grid data; Step S353: Calculate the local stress gradient data based on the initial transition region grid data and the stress-strain relationship data; Step S354: Perform adaptive optimization on the transition region grid in the initial transition region grid data based on the local stress gradient data. When the stress gradient is greater than the preset threshold, subdivide the grid to 0.5 mm to obtain the optimized grid data of the transition region; Step S355: Calculate the stress distribution in the transition region based on the optimized grid data and the stress-strain relationship data, so as to obtain the third deformation field data of the transition region.

8. The method for three-dimensional image reconstruction and positioning of lung puncture according to claim 7, wherein, Step S4 includes the following steps: Step S41: Calculate the gray gradient of adjacent images in the CT scan sequence data to obtain the image gray gradient data, where the image gray gradient data includes the x-direction gradient data, the y-direction gradient data, and the time-direction gradient data; Step S42: Construct the initial motion vector field data based on the first deformation field data of the central region, the second deformation field data of the peripheral region, and the third deformation field data of the transition region; Step S43: Establish a brightness constancy constraint equation according to the image gray gradient data to obtain the brightness constraint condition data; Step S44: Establish a velocity smoothness constraint equation according to the initial motion vector field data to obtain the velocity constraint condition data; Step S45: Perform weighted combination on the brightness constraint condition data and the velocity constraint condition data to obtain the optical flow method energy functional expression; and perform Euler-Lagrange equation derivation on the optical flow method energy functional expression to obtain the control equation set of the variational problem; Step S46: Numerically solve the control equation set to obtain the optimized motion vector field data of each region; Step S47: Perform data fusion according to the optimized motion vector field data of each region, eliminate the discontinuity between regions through boundary smoothing processing, construct an error function and compare it with the preset threshold. When the error is greater than the threshold, return to Step S45 to adjust the weight coefficient and recalculate until the error is less than the threshold, so as to obtain the tissue deformation field data.

9. The method for three-dimensional image reconstruction and positioning of lung puncture according to claim 8, wherein, Step S5 includes the following steps: Step S51: Perform grayscale threshold segmentation on the CT sequence data to obtain initial lung contour data; extract lung surface feature line data based on the initial lung contour data, where the lung surface feature line data includes the surface projections of the pulmonary fissures and the main vascular courses. Step S52: Perform temporal analysis on the tissue deformation field data to obtain tissue deformation law data for different respiratory phases; establish an elastic deformation model based on the tissue deformation law data using the bubble animation simulation algorithm, regarding the tissue as a droplet with surface tension, so as to obtain surface tension parameter data. Step S53: Construct a dynamic equation describing lung tissue deformation based on the surface tension parameter data and the tissue deformation law data, so as to obtain deformation dynamics model data. Step S54: Numerically solve the deformation dynamics model data to obtain three-dimensional grid node position data for different respiratory phases; perform surface reconstruction based on the three-dimensional grid node position data and maintain the original grid density in the feature line region, so as to obtain initial three-dimensional model data. Step S55: Simplify the mesh of the initial three-dimensional model data in the non-feature line region, reducing the number of meshes to 50% of the original through the quadratic error metric criterion, so as to obtain optimized three-dimensional model data; identify the spatial distributions of the ribs and major blood vessels in the optimized three-dimensional model data, so as to obtain forbidden penetration region data. Step S56: Calculate feasible puncture path sets based on the forbidden penetration region data and the respiratory phase mapping relationship, so as to obtain initial path set data; perform safety assessment on the initial path set data to obtain path risk assessment data, where the safety assessment specifically calculates the minimum distance between each path and the forbidden penetration region for different respiratory phases. Step S57: Screen and sort the initial path set according to the path risk assessment data to obtain the final puncture path planning model data.

10. A three-dimensional image reconstruction and positioning system for lung puncture, characterized in that, A system for three-dimensional image reconstruction and localization for lung puncture, which is used to execute the method for three-dimensional image reconstruction and localization for lung puncture as described in Claim 1, and the system for three-dimensional image reconstruction and localization for lung puncture includes: A data acquisition module, which is used to acquire the CT scan sequence data and multi-point chest wall infrared marker data of the patient, where the multi-point chest wall infrared marker data includes respiratory reference position data, real-time angular velocity data of the chest wall movement state, and real-time acceleration data of each point on the chest wall. A tissue mechanics modeling module, which is used to predict the non-linear viscoelasticity of lung tissue based on the multi-point chest wall infrared marker data using a Maxwell model, and input the real-time angular velocity data and the real-time acceleration data into the non-linear viscoelasticity prediction model of lung tissue for iterative calculation to obtain lung tissue stress-strain relationship data, where the stress-strain relationship data includes instantaneous stress response data and delayed strain component data. A partitioned mesh generation module, which is used to perform partitioned mesh generation on the lung tissue based on the stress-strain relationship data to obtain first deformation field data of the central region based on high-density meshes, second deformation field data of the peripheral region based on low-density meshes, and third deformation field data of the transition region constructed using adaptive meshes based on local stress gradients in the remaining region. A motion vector calculation module, which is used to calculate the motion vector field for the first deformation field data in the central region, the second deformation field data in the peripheral region, and the third deformation field data in the transition region by using the optical flow method, construct an energy functional according to the brightness constancy constraint and the velocity smoothness constraint, and solve it by the variational method to obtain the optimized tissue deformation field data; A three-dimensional reconstruction and path planning module, which is used to apply the tissue deformation field data to the CT scan sequence data and perform three-dimensional reconstruction based on the bubble animation simulation algorithm, simplify the grid in the non-feature line region while maintaining the grid density in the feature line region, and obtain the puncture path planning model data with respiratory phase mapping.

Images