Real-time positioning and puncture trajectory control method and system in breast ductoscopy surgery

Abstract

The application provides a real-time positioning and puncture trajectory control method and system in a breast duct endoscopy operation, relates to the technical field of medical devices, and comprises the following steps: establishing a real-time coordinate system by collecting endoscopy image sequences and electromagnetic position data, obtaining three-dimensional depth information of a lumen by using a light field reconstruction algorithm, adjusting and calculating a stress field distribution in a grid form, determining an optimal trajectory based on a puncture safety coefficient, and adjusting and compensating parameters in real time by using a fiber strain sensing array. The application realizes accurate positioning and trajectory control of a puncture device in a breast duct endoscopy operation, and improves the safety and accuracy of the operation.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of medical devices, in particular to a real-time positioning and puncture trajectory control method and system in a breast ductoscopy operation. BACKGROUND

[0002] Breast ductoscopy is a minimally invasive technique for examining the internal conditions of breast ducts, widely used in early diagnosis and treatment of breast diseases. In breast ductoscopy surgery, doctors need to enter the internal ducts of the breast through a small endoscope probe to observe, biopsy or treat suspicious lesion tissues. With the development of medical technology, breast ductoscopy technology has developed from a simple diagnostic tool to a comprehensive platform that can perform biopsy and minimally invasive treatment.

[0003] The current positioning and puncture technology in breast ductoscopy surgery faces many challenges and deficiencies. The traditional breast ductoscopy positioning method mainly relies on the experience of doctors and two-dimensional image information for judgment, lacking precise three-dimensional spatial positioning ability, making it difficult for doctors to accurately locate the lesion site, especially in areas with complex breast duct branches, leading to positioning deviation and missed diagnosis.

[0004] The existing breast ductoscopy puncture control technology fails to effectively solve the tissue deformation problem. When the puncture instrument contacts the breast duct tissue, the soft tissue will deform, causing deviation between the pre-planned puncture trajectory and the actual operation, reducing the accuracy of puncture and increasing the risk of damage to surrounding healthy tissue. The existing technology lacks a real-time monitoring mechanism for the force of the instrument during puncture, and doctors cannot obtain accurate force feedback information from the puncture needle and the tissue, making it difficult to judge the tissue boundary and puncture depth during puncture, which may cause excessive puncture or insufficient puncture, affecting the safety and effectiveness of the operation. SUMMARY

[0005] The present application provides a real-time positioning and puncture trajectory control method and system in a breast ductoscopy operation, which can solve the problems in the prior art.

[0006] In a first aspect, the present application provides a real-time positioning and puncture trajectory control method in a breast ductoscopy operation, comprising:

[0007] Collecting an endoscopic image sequence of a breast ductoscopy probe and obtaining spatial positioning data through an electromagnetic position sensor on the probe, time-series registering the endoscopic image sequence and the spatial positioning data, and establishing a real-time coordinate system of the operation space;

[0008] Based on the real-time coordinate system, processing the endoscopic image sequence through a light field reconstruction algorithm to obtain three-dimensional depth information of the inner wall of the lumen, and combining the spatial positioning data of the electromagnetic position sensor to determine the real-time deformation parameters of the lumen channel;

[0009] Adopting a grid form adjustment method to deform the lumen wall surface adaptively, using strain energy gradient to control node optimization, calculating the stress field distribution of the lumen wall surface according to the real-time deformation data;

[0010] According to the stress field distribution, calculating the puncture safety coefficient through multi-layer sampling points, determining the optimal puncture trajectory based on the puncture safety coefficient, and realizing the motion tracking control of the puncture instrument combined with the speed control point sequence;

[0011] Using an optical fiber strain sensing array to collect the deformation signal of the puncture instrument in real time, and calculating the real-time force data of the puncture instrument;

[0012] According to the real-time force data, analyzing the deviation between the actual motion trajectory of the puncture instrument and the puncture trajectory, adjusting the compensation parameters in real time, and realizing the precise trajectory tracking control of the puncture instrument.

[0013] In an alternative embodiment, based on the real-time coordinate system, the endoscopic image sequence is processed by a light field reconstruction algorithm to obtain the three-dimensional depth information of the lumen inner wall, and the real-time deformation parameters of the lumen channel are determined combined with the spatial positioning data of the electromagnetic position sensor, including:

[0014] Setting up a light field sensor array in the real-time coordinate system, collecting multi-view light information of the lumen inner wall through the light field sensor array, and obtaining the principal ray coordinates and micro ray coordinates by light field reconstruction processing on the multi-view light information;

[0015] According to the principal ray coordinates and the micro ray coordinates, calculating the light space re-projection error, constructing an optimization objective combined with a depth smoothing constraint term, and obtaining the three-dimensional depth information of the lumen inner wall by iterative optimization through gradient descent method;

[0016] Performing coordinate transformation on the spatial positioning data of the electromagnetic position sensor and the three-dimensional depth information to generate three-dimensional point cloud data representing the lumen inner wall;

[0017] Performing feature matching and registration on the continuously collected multi-frame three-dimensional point cloud data, extracting point cloud deformation vectors, and obtaining the real-time deformation parameters of the lumen channel by introducing elastic constraints and iterative optimization.

[0018] In an alternative embodiment, performing feature matching and registration on the continuously collected multi-frame three-dimensional point cloud data, extracting point cloud deformation vectors, and obtaining the real-time deformation parameters of the lumen channel by introducing elastic constraints and iterative optimization, including:

[0019] Calculating the point cloud curvature, normal vector and geometric moment distribution of the continuously collected multi-frame three-dimensional point cloud data to obtain feature descriptors representing the local geometric characteristics of each point cloud;

[0020] Based on the feature descriptor, a K nearest neighbor search algorithm is used to perform feature matching on the adjacent frame point cloud to obtain initial matching point pairs, and a transformation relationship between the initial matching point pairs is calculated through iteration, combined with abnormal matching point pair elimination, to obtain a point cloud registration result;

[0021] Based on the point cloud registration result, a deformation direction and a deformation amplitude of a local point cloud are extracted, a point cloud deformation vector is constructed, and a deformation distribution feature is obtained through gradient analysis of the point cloud deformation vector;

[0022] Based on the point cloud deformation vector transmission of the nearest neighbor point pair, a regional deformation constraint is constructed, and a corresponding elastic correction value is calculated. According to the elastic correction value, a deformation gradient and an elastic deformation variable of the adjacent region of the point cloud are extracted, a moving least squares method is used to establish a deformation parameter mapping model, and the point cloud deformation vector is iteratively optimized based on the deformation parameter mapping model. The real-time deformation parameter of the lumen channel is determined by calculating the difference value of the point cloud deformation vector of adjacent iteration rounds.

[0023] In an optional embodiment, a grid morphology adjustment method is used to perform adaptive deformation on the lumen wall surface, a strain energy gradient is used to control node optimization, and a stress field distribution of the lumen wall surface is calculated according to the real-time deformation data, including:

[0024] Strain field information in the real-time deformation data is extracted, and a mapping relationship between the strain field information and the nodes of the grid element is established;

[0025] The strain rate tensor of the grid element is calculated according to the strain field information, and the principal strain direction and the principal strain rate of the grid element are obtained through eigenvalue decomposition of the strain rate tensor;

[0026] Deformation control parameters of the grid element are constructed based on the principal strain direction and the principal strain rate; and the node coordinates of the grid element are adjusted according to the deformation control parameters;

[0027] The strain energy gradient between the adjusted adjacent grid elements is calculated, and the node position iterative optimization is triggered when the strain energy gradient exceeds a preset gradient threshold, until the strain energy gradient is less than or equal to the gradient threshold;

[0028] A stress interpolation function is established according to the optimized grid node distribution, the order of the stress interpolation function is positively correlated with the strain rate of the corresponding grid element, and the stress field distribution of the lumen wall surface is calculated.

[0029] In an optional embodiment, constructing deformation control parameters of the grid element based on the principal strain direction and the principal strain rate includes:

[0030] The deformation control parameters include a tensile coefficient along the principal strain direction and a compression coefficient perpendicular to the principal strain direction; a ratio of a principal strain rate of the grid cell to a preset reference strain rate is calculated, the ratio is taken as a tensile coefficient reference value, the tensile coefficient along the principal strain direction is determined according to the tensile coefficient reference value, and the tensile coefficient non-linearly increases with an increase of the ratio;

[0031] A deviation of an aspect ratio of the grid cell from a preset shape factor is calculated, and the compression coefficient perpendicular to the principal strain direction is determined based on the deviation, and the compression coefficient is used to constrain grid distortion;

[0032] The node coordinates of the grid cell are adjusted according to the deformation control parameters, including:

[0033] The node coordinates of the grid cell are projected to a local coordinate system with the principal strain direction as a reference, in the local coordinate system, a displacement amount of a node along the principal strain direction is proportional to the tensile coefficient, and a displacement amount of the node perpendicular to the principal strain direction is proportional to the compression coefficient;

[0034] The node coordinates of the grid cell in the global coordinate system are updated based on the displacement amounts.

[0035] In an optional embodiment, a puncture safety coefficient is calculated through the multi-layer sampling points according to the stress field distribution, an optimal puncture trajectory is determined based on the puncture safety coefficient, and the motion tracking control of the puncture instrument is realized in combination with a speed control point sequence, including:

[0036] Multi-layer sampling points are constructed in the surgical space according to the stress field distribution, and each layer of sampling points is connected to form a sampling plane; a connection relationship is established between adjacent sampling planes, a puncture safety coefficient between the sampling points is calculated according to the connection relationship, and the puncture safety coefficient is inversely proportional to a stress field value at the sampling point;

[0037] A starting point and a target point of puncture are set, a virtual potential field distribution map based on the puncture safety coefficient is constructed, a resultant force direction of an attractive force and a repulsive force is calculated in the virtual potential field distribution map, and an optimal puncture path is calculated through adaptive step gradient iteration;

[0038] Speed control points are set along the optimal puncture path, and a motion speed of the puncture instrument at each speed control point is determined according to a change in curvature of the lumen;

[0039] A tangent vector and a normal vector are obtained at each speed control point, and a motion posture parameter of the puncture instrument is calculated;

[0040] Real-time motion parameters of the puncture instrument are collected, the real-time motion parameters are compared with the motion posture parameter, and a deviation compensation signal is generated;

[0041] Adjust the driving output of the puncture instrument according to the deviation compensation signal, so as to realize continuous motion control of the puncture instrument along the optimal puncture path.

[0042] In an optional embodiment, the starting point and target point of puncture are set, a virtual potential field distribution map based on the puncture safety coefficient is constructed, the direction of the resultant force of the attractive force and the repulsive force is calculated in the virtual potential field distribution map, and the optimal puncture path is calculated through adaptive step gradient iteration, comprising:

[0043] The spatial coordinates of the puncture starting point and target point are determined on the sampling plane, a basic potential field function is constructed with the reciprocal of the puncture safety coefficient, a local potential field distribution is established at each sampling point, and the local potential fields of adjacent sampling planes are interpolated and fitted to form a continuous virtual potential field distribution map;

[0044] An attractive force field is established with the target point as the attractive force source, the strength of the attractive force field decays with the increase of the distance, a repulsive force field is established with the sampling points whose puncture safety coefficients are less than a preset safety threshold as the repulsive force source, the strength of the repulsive force field is inversely proportional to the puncture safety coefficient, the vector sum of the attractive force and the repulsive force is superimposed at the calculation points of the virtual potential field distribution map, and the direction of the resultant force is determined;

[0045] An adaptive step gradient descent algorithm is adopted, the starting point is taken as the initial search position, iterative calculation is performed along the direction of the resultant force until the position error between the search position and the target point is less than a preset error threshold, the optimal puncture path is obtained, and the optimal puncture path is smoothed to ensure the continuity of the curvature.

[0046] In a second aspect of the embodiment of the present application, a real-time positioning and puncture trajectory control system in a breast ductoscopy surgery is provided, comprising:

[0047] A first unit is configured to collect an endoscopic image sequence of a breast ductoscopy probe, acquire spatial positioning data through an electromagnetic position sensor on the probe, time sequence register the endoscopic image sequence and the spatial positioning data, and establish a real-time coordinate system of a surgery space;

[0048] A second unit is configured to process the endoscopic image sequence through a light field reconstruction algorithm based on the real-time coordinate system, acquire three-dimensional depth information of a lumen inner wall, determine real-time deformation parameters of a lumen channel in combination with the spatial positioning data of the electromagnetic position sensor, and calculate a stress field distribution of the lumen wall surface according to the real-time deformation data.

[0049] A third unit is configured to perform adaptive deformation on the lumen wall surface by using a grid shape adjustment method, control node optimization by using a strain energy gradient, and calculate the stress field distribution of the lumen wall surface according to the real-time deformation data.

[0050] The fourth unit is configured to calculate a puncture safety coefficient through a plurality of sampling points according to the stress field distribution, determine an optimal puncture trajectory based on the puncture safety coefficient, and realize motion tracking control of the puncture instrument in combination with a speed control point sequence.

[0051] The fifth unit is configured to collect deformation signals of the puncture instrument in real time by using a fiber strain sensing array, and calculate real-time stress data of the puncture instrument.

[0052] The sixth unit is configured to analyze deviation between an actual motion trajectory of the puncture instrument and the puncture trajectory according to the real-time stress data, adjust compensation parameters in real time, and realize precise trajectory tracking control of the puncture instrument.

[0053] In a third aspect, an electronic device is provided, including:

[0054] a processor;

[0055] a memory for storing processor-executable instructions;

[0056] The processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0057] In a fourth aspect, a computer-readable storage medium is provided, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method described above.

[0058] In the embodiment of the present application, the real-time positioning and puncture trajectory control method in the breast ductoscopy surgery is provided, fusion processing of the breast ductoscopy endoscopy image and the electromagnetic position sensing data is collected, a dynamic real-time coordinate system is established, precise three-dimensional reconstruction of the inner wall of the breast duct cavity is realized, the surgical navigation positioning is more accurate and reliable, and the problem of inaccurate positioning in the traditional breast ductoscopy surgery is effectively solved; the duct cavity wall surface modeling method based on the light field reconstruction and the strain energy gradient is adopted, real-time deformation monitoring and stress field analysis are combined, a puncture safety coefficient evaluation mechanism is constructed, the optimal puncture path can be scientifically and reasonably planned, the safety of the puncture operation is greatly improved, the risk of injury is reduced, and the method is particularly suitable for precise intervention of complex lesion sites; the stress state of the puncture instrument is monitored in real time by the fiber strain sensing array, a closed-loop feedback regulation mechanism is established, the puncture trajectory can be dynamically adjusted during the surgery, deviation caused by tissue deformation is compensated, the puncture positioning accuracy is significantly improved, the surgical effect is improved, complications are reduced, and the method has important value for improving the treatment effect of breast minimally invasive surgery. BRIEF DESCRIPTION OF DRAWINGS

[0059] Figure 1 FIG. 1 is a flowchart of the real-time positioning and puncture trajectory control method in the breast ductoscopy surgery according to the embodiment of the present application;

[0060] Figure 2 a heat map for point cloud deformation vector analysis;

[0061] Figure 3 an error scatter plot for lumen wall surface stress field calculation;

[0062] Figure 4 a rendering for virtual potential field-based puncture path planning. DETAILED DESCRIPTION

[0063] To make the objects, technical solutions and advantages of the embodiments of the present application clearer, the following will be combined with the drawings in the embodiments of the present application to make a clear and complete description of the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0064] The technical solutions of the present application will be described in detail in the following specific embodiments. The following specific embodiments can be combined with each other, and the same or similar concepts or processes can not be described in some embodiments.

[0065] Figure 1 A flowchart of a real-time positioning and puncture trajectory control method in a mammary ductoscopy operation according to an embodiment of the present application is shown in FIG. 1, which includes the following steps: Figure 1

[0066] Collecting an endoscopic image sequence of a mammary ductoscopy probe and acquiring spatial positioning data through an electromagnetic position sensor on the probe, time-series registering the endoscopic image sequence and the spatial positioning data, and establishing a real-time coordinate system of the operation space;

[0067] Based on the real-time coordinate system, processing the endoscopic image sequence through a light field reconstruction algorithm to acquire three-dimensional depth information of the inner wall of the lumen, and combining the spatial positioning data of the electromagnetic position sensor to determine real-time deformation parameters of the lumen channel;

[0068] Adopting a grid shape adjustment method to adaptively deform the lumen wall surface, using a strain energy gradient to control node optimization, and calculating stress field distribution of the lumen wall surface according to the real-time deformation data;

[0069] According to the stress field distribution, calculating a puncture safety coefficient through multi-layer sampling points, determining an optimal puncture trajectory based on the puncture safety coefficient, and realizing movement tracking control of a puncture instrument in combination with a speed control point sequence;

[0070] Using a fiber strain sensing array to collect deformation signals of the puncture instrument in real time, and calculating real-time force data of the puncture instrument;

[0071] ​According to the real-time force data, a deviation between an actual motion trajectory of the puncture instrument and a puncture trajectory is analyzed, a compensation parameter is adjusted in real time, and precise trajectory tracking control of the puncture instrument is realized.

[0072] In an optional embodiment, based on the real-time coordinate system, an endoscopic image sequence is processed by a light field reconstruction algorithm to obtain three-dimensional depth information of the intraluminal wall, and real-time deformation parameters of the lumen passage are determined in combination with spatial positioning data of the electromagnetic position sensor, including:

[0073] A light field sensor array is arranged in the real-time coordinate system, multi-view light information of the intraluminal wall is collected by the light field sensor array, and main light coordinates and micro light coordinates are obtained by light field reconstruction processing on the multi-view light information;

[0074] Light space re-projection errors are calculated according to the main light coordinates and the micro light coordinates, an optimization target is constructed in combination with a depth smoothing constraint term, and three-dimensional depth information of the intraluminal wall is obtained by iterative optimization of the gradient descent method;

[0075] The spatial positioning data of the electromagnetic position sensor and the three-dimensional depth information are subjected to coordinate transformation to generate three-dimensional point cloud data representing the intraluminal wall;

[0076] Feature matching and registration are performed on the continuously collected multiple frames of three-dimensional point cloud data, point cloud deformation vectors are extracted, and real-time deformation parameters of the lumen passage are obtained by introducing elastic constraints and iterative optimization.

[0077] In a specific embodiment, a light field sensor array is arranged in a real-time coordinate system, the light field sensor array includes a plurality of microlens units, and each microlens unit is provided below with a plurality of photosensitive elements. The light field sensor array is inserted into the lumen by an endoscope to collect multi-view light information of the intraluminal wall. Specifically, when light is projected onto the photosensitive element array through the microlens array, the photosensitive element under each microlens records light information incident at different angles. Assuming that the light field sensor array used contains 11x11 microlens units, and each microlens is provided below with 7x7 photosensitive elements, 121 different view image information can be obtained.

[0078] The collected multi-view light information is processed by light field reconstruction to obtain the main light ray coordinates and the micro light ray coordinates. The main light ray coordinates represent the position of the light ray passing through the center of the micro lens, usually represented by (s, t); the micro light ray coordinates represent the projection position of the light ray on the image plane, usually represented by (u, v). In actual processing, for the light intensity value I recorded by each photosensitive element, it can be mapped into the four-dimensional light field function L(s, t, u, v). For example, when the micro lens center position is (3, 4) and the light ray projection position is (2, 1), it can be represented as L(3, 4, 2, 1) = I, where I is the light intensity value recorded at this position, and the value range is 0-255.

[0079] The light ray space reprojection error is calculated according to the main light ray coordinates and the micro light ray coordinates. Specifically, for each point P on the intraluminal wall, the depth value of which is d, the projection position of the point at different viewing angles can be calculated through the light field imaging model. If the depth value d is accurate, the projection positions at all viewing angles should be consistent. By comparing the difference between the actually observed light ray and the theoretical projection position, the reprojection error E is calculated. In actual application, for a point P on the intraluminal wall, assuming that the initial depth estimation value of the point is 5 mm, it is found through calculation that the average reprojection error at 11 different viewing angles is 2.3 pixels, indicating that the depth estimation is not accurate enough and needs to be optimized.

[0080] An optimization target is constructed in combination with a depth smoothing constraint term. The depth smoothing constraint is based on the feature that the intraluminal wall surface usually has local continuity, and requires that the depth values of adjacent points should change smoothly. The optimization target function J includes the reprojection error term E and the depth smoothing constraint term S, and the two are balanced by a weight factor λ (taking a value of 0.5). The depth value d is iteratively optimized by the gradient descent method to minimize the target function J. In the actual optimization process, the learning rate α is set to 0.01, and the maximum number of iterations is set to 200 times. When the target function value changes less than 0.001 for 5 consecutive iterations, it is considered that the optimization converges. After optimization, the final depth value of the above example point P is optimized to 6.2 mm, and the reprojection error is reduced to 0.8 pixels.

[0081] The spatial positioning data of the electromagnetic position sensor is coordinate-transformed with the three-dimensional depth information to generate three-dimensional point cloud data representing the inner wall of the lumen. The electromagnetic position sensor is installed at the tip of the endoscope and can provide position and attitude information of the endoscope in an external reference coordinate system, including position coordinates (X, Y, Z) and a rotation matrix R. Through a known rigid transformation relationship T between the endoscope coordinate system and the light field sensor coordinate system, combined with the depth information d obtained by light field reconstruction, the point in the light field coordinate system can be converted into the external reference coordinate system. For example, the position of the tip of the endoscope at a certain time is (120 mm, 45 mm, 80 mm), and the attitude matrix is R (this 3x3 matrix describes the orientation of the device), combined with the depth information, the light field coordinates (3, 4, 6.2 mm) of the point P on the inner wall of the lumen can be converted into a three-dimensional point (126.5 mm, 48.2 mm, 85.7 mm) in the global coordinate system.

[0082] Feature matching and registration are performed on the continuously collected multiple frames of three-dimensional point cloud data. Feature points such as corner points and edge points are extracted from the surface of the inner wall of the lumen, and a corresponding relationship between the feature points of the front and rear frames is established. The iterative closest point (ICP) algorithm is used for point cloud registration to align the point clouds between the consecutive frames. In actual application, to improve the registration accuracy, about 200 feature points can be extracted from each frame of point cloud, and the corresponding relationship is established through feature descriptor matching. For example, in the point cloud data collected at two consecutive times, the feature point matching success rate reaches 85%, and the root mean square error after registration is 0.3 mm.

[0083] The point cloud deformation vector is extracted, and the real-time deformation parameters of the lumen passage are obtained by introducing elastic constraints and iterative optimization. The point cloud deformation vector is defined as the displacement vector between the corresponding points after registration. Considering the elastic properties of biological tissues, elastic constraints are introduced to limit the reasonable range of deformation. Through an iterative optimization algorithm, the real-time deformation parameters of the lumen passage are calculated, including local stretch rate, compression rate, and bending degree. In a certain test, the maximum radial deformation of the lumen under respiratory motion was detected to be 2.1 mm, the axial deformation was 3.8 mm, and the local maximum stretch rate was 8.5%. These parameters can be used to guide the design of subsequent medical devices and surgical planning.

[0084] In an alternative embodiment, feature matching and registration are performed on the continuously collected multiple frames of three-dimensional point cloud data, the point cloud deformation vector is extracted, and the real-time deformation parameters of the lumen passage are obtained by introducing elastic constraints and iterative optimization, including:

[0085] The point cloud curvature, normal vector, and geometric moment distribution of the continuously collected multiple frames of three-dimensional point cloud data are calculated to obtain feature descriptors representing the local geometric characteristics of each point cloud;

[0086] Based on the feature descriptor, K-Nearest Neighbor search algorithm is used to perform feature matching on the adjacent frame point cloud to obtain initial matching point pairs, and through iterative calculation of the transformation relationship between the initial matching point pairs, combined with the elimination of abnormal matching point pairs, a point cloud registration result is obtained.

[0087] Based on the point cloud registration result, the deformation direction and deformation amplitude of the local point cloud are extracted, the point cloud deformation vector is constructed, and the deformation distribution characteristics are obtained through gradient analysis of the point cloud deformation vector.

[0088] Based on the point cloud deformation vector transmission of the nearest neighbor point pair, a regional deformation constraint is constructed, and a corresponding elastic correction value is calculated. According to the elastic correction value, the deformation gradient and elastic deformation variable of the adjacent region of the point cloud are extracted, a moving least squares method is used to establish a deformation parameter mapping model, and based on the deformation parameter mapping model, the point cloud deformation vector is iteratively optimized for multiple rounds. The real-time deformation parameters of the lumen channel are determined by calculating the difference value of the point cloud deformation vectors of adjacent iteration rounds.

[0089] In a specific embodiment, in the point cloud feature descriptor calculation stage, for each frame of three-dimensional point cloud data, the curvature information of the point cloud is calculated. Specifically, for each point P in the point cloud, all points in a spherical neighborhood with a radius of 3mm are selected to form a point set N(P). By performing principal component analysis on the points in N(P), three eigenvalues of the covariance matrix λ1≥λ2≥λ3≥0 can be obtained. The curvature of point P can be represented by the formula λ3 / (λ1+λ2+λ3). In practical applications, when the value is less than 0.02, it is considered that the point is located in a planar region; when the value is greater than 0.05, it is considered that the point is located in an edge or corner region. The normal vector of the point, i.e. the eigenvector corresponding to the smallest eigenvalue of the covariance matrix, is calculated. In order to enhance the robustness of the feature description, the geometric moment distribution of the point cloud is further calculated, including dividing the space into 8 octants in the neighborhood of point P, and calculating the proportion of the number of points in each octant to form an 8-dimensional feature vector. Combining the curvature value, normal vector and geometric moment distribution, a 13-dimensional feature descriptor is constructed to represent the local geometric characteristics of the point cloud.

[0090] In the point cloud feature matching and registration stage, for two consecutive frames of point cloud data Frame_i and Frame_i+1, K-Nearest Neighbor search algorithm is used for feature matching. For each point in Frame_i, K points (K=5) with the smallest Euclidean distance of feature descriptors in Frame_i+1 are searched, and the point pairs with a distance less than a preset threshold (1.5 times the average distance of feature descriptors) are taken as initial matching point pairs. In order to eliminate false matches, a Random Sample Consensus algorithm is used for optimization. Randomly select 3 pairs of points from the initial matching point pairs, calculate the rigid transformation matrix, and then apply the transformation matrix to all points in Frame_i, calculate the distance between the transformed points and the corresponding points in Frame_i+1. If the distance is less than 1 mm, it is considered that the point pair is consistent. Repeat this process 1000 times, and select the transformation matrix with the most consistent point pairs as the optimal solution. Based on the optimal transformation matrix, the two frames of point cloud are preliminarily registered, and then the Iterative Closest Point algorithm is used for fine registration, and finally the point cloud registration result is obtained, and the registration accuracy can reach within 0.5 mm.

[0091] In the point cloud deformation vector extraction stage, based on the registration result, for each point P in the registered Frame_i, find the nearest neighbor point Q in Frame_i+1. Calculate the vector PQ as the deformation vector v(P) of point P. In order to eliminate the influence of noise, the deformation vector is smoothed, and the deformation vector of each point is weighted and averaged with the deformation vectors of the points in its neighborhood (radius 2 mm), and the weight is inversely proportional to the distance. Gradient analysis is performed on the point cloud deformation vector, the spatial derivative of the deformation vector field is calculated, and the deformation gradient tensor is obtained. In practical application, when the direction consistency of the deformation vector in the local region is higher than 0.85 (cosine similarity), and the amplitude change is less than 20%, it is determined that the region is a rigid deformation region; when the direction consistency is between 0.5 and 0.85, and the amplitude change is between 20% and 50%, it is determined as an elastic deformation region; otherwise, it is determined as a complex deformation region.

[0092] In the elastic constraint and deformation parameter iterative optimization stage, the regional deformation constraint is constructed based on the point cloud deformation vector. For each point P in the point cloud, its influence domain is defined as all points within a radius r (r = 5mm). For each point Q in the influence domain, its deformation constraint is that v(Q) should have a similar direction and gradually changing amplitude as v(P). Specifically, the elastic correction value δv(Q) is calculated so that the corrected deformation vector v'(Q) = v(Q) + δv(Q) satisfies the regional deformation constraint. The calculation of the elastic correction value takes into account the distance factor between points, the closer the distance, the greater the constraint force. Based on the elastic correction value, the deformation gradient and elastic deformation variable of the adjacent region are further extracted. The moving least squares method is used to establish a deformation parameter mapping model, which can smoothly map the deformation parameters of discrete points to the entire space. In the model construction, feature points (such as points with large curvature changes) are selected as control points, and each control point has an influence range of 10mm radius. Based on the deformation parameter mapping model, the point cloud deformation vector is iteratively optimized for multiple rounds. In each iteration, the deformation vector of the point is updated to better meet the regional deformation constraint. When the average difference value of the deformation vector between two adjacent iterations is less than 0.1mm, it is considered that the optimization converges, and the final real-time deformation parameter of the lumen channel is obtained.

[0093] In the representation of the deformation parameter, the lumen channel is divided into several cross sections, and 16 points are uniformly sampled on each cross section. The radial displacement and tangential displacement of these points are calculated to represent the expansion / contraction and bending / twisting of the lumen.

[0094] The prior art has obvious deficiencies in real-time positioning and puncture trajectory control in the breast ductoscopy surgery. Traditional methods usually use point cloud registration algorithms based on rigid transformation, such as the Iterative Closest Point (ICP) algorithm. This kind of method assumes that there is only a rigid transformation between two frames of point clouds, and cannot effectively handle the elastic deformation of breast tissue during the surgery. In terms of feature extraction, the prior art mostly uses simple geometric features (such as curvature, normal vector) or histogram-based descriptors (such as FPFH, SHOT, etc.). These descriptors have low matching accuracy when dealing with low-texture and high-similarity biological tissues, which can easily lead to registration failure. In terms of point cloud deformation parameter optimization, the prior art mostly uses global optimization methods, which have high computational complexity and are difficult to meet the real-time requirements of the surgery, and lack of constraint mechanisms for the deformation characteristics of biological tissues.

[0095] As Figure 2The point cloud local deformation direction and deformation amplitude distribution extracted by the method of the embodiment are shown. By constructing the point cloud deformation vector and performing gradient analysis, the heat map clearly presents the deformation differences in different regions. The central red area in the figure represents the maximum deformation area, with a maximum deformation of 5.2 mm. The center point as the origin radiates outward, showing a clear gradient distribution. The deformation vector arrow direction and length around the center intuitively represent the deformation direction and amplitude. From the heat map, it can be seen that the deformation amplitude gradually decreases from the center to the outside, and in the edge region, it is only 1.2-1.6 mm. By analyzing the deformation gradient in the heat map, it can be calculated that the deformation gradient range is 0.3-1.2 mm / cm, and the total deformation area is 76.2 cm 2 The detection sensitivity of the method is 0.2 mm, which is much better than the traditional method. Based on the deformation heat map, the system calculates the average elastic correction value as 0.38 mm, mainly concentrated in the range of 3.5 cm around the center point. This accurate deformation analysis benefits from the deformation parameter mapping model established by the method of the embodiment using the moving least squares method, which realizes high-precision deformation quantization through 8 rounds of iterative optimization. Compared with the traditional point cloud analysis method, the method of the embodiment not only can capture the tiny deformation, but also can accurately represent the directionality and regional distribution of the deformation, providing accurate parameter basis for subsequent elastic correction. It has important significance for real-time deformation monitoring of structures such as lumen channels.

[0096] The method of the embodiment proposes a real-time positioning and puncture trajectory control method suitable for breast ductoscopy surgery. In the feature descriptor design, the point cloud curvature, normal vector and geometric moment distribution are combined to construct a 13-dimensional feature descriptor, which significantly improves the recognition ability of feature points on the biological tissue surface. In the point cloud registration, unlike the traditional rigid registration algorithm, an abnormal matching point elimination mechanism based on random sample consensus is introduced, which effectively improves the registration accuracy, so that the registration accuracy is within 0.5 mm, meeting the accuracy requirements of breast ductoscopy surgery. In the deformation analysis, a deformation vector optimization method based on local area elastic constraint is proposed. By defining the deformation constraint relationship of points in the influence domain, the elastic correction value is calculated, effectively capturing the non-rigid deformation characteristics of breast tissue. In the deformation parameter optimization, the moving least squares method is used to establish the deformation parameter mapping model, and through multiple rounds of iterative optimization, the smooth estimation of the deformation parameter is realized, while ensuring the calculation efficiency, supporting real-time analysis and control during the surgery. Under the condition of tissue deformation, the positioning accuracy is maintained, and the deformation parameters required for puncture trajectory control are obtained in real time. By introducing elastic constraint and deformation gradient analysis, the method of the embodiment can accurately distinguish between rigid deformation area, elastic deformation area and complex deformation area, providing a theoretical basis for accurate control of the puncture trajectory.

[0097] In an alternative embodiment, the lumen wall surface is adaptively deformed using a mesh morphing method, and a strain energy gradient control node optimization is used to calculate the stress field distribution of the lumen wall surface according to the real-time deformation data, which includes:

[0098] Strain field information is extracted from the real-time deformation data, and a mapping relationship between the strain field information and the nodes of the mesh elements is established.

[0099] The strain rate tensor of the mesh element is calculated according to the strain field information, and the principal strain direction and the principal strain rate of the mesh element are obtained by eigenvalue decomposition of the strain rate tensor.

[0100] Deformation control parameters of the mesh element are constructed based on the principal strain direction and the principal strain rate, and the node coordinates of the mesh element are adjusted according to the deformation control parameters.

[0101] The strain energy gradient between the adjusted adjacent mesh elements is calculated, and the node position iterative optimization is triggered when the strain energy gradient exceeds a preset gradient threshold, until the strain energy gradient is less than or equal to the gradient threshold.

[0102] A stress interpolation function is established according to the optimized mesh node distribution, the order of the stress interpolation function is positively correlated with the strain rate of the corresponding mesh element, and the stress field distribution of the lumen wall surface is calculated.

[0103] In a specific embodiment, real-time deformation data of the lumen wall surface is obtained, which can be obtained by medical imaging devices such as ultrasound, CT or MRI. In a specific case, a blood vessel is modeled, and the obtained lumen deformation data contains spatial displacement information of 500 discrete points, each point has x, y, z three-dimensional displacement values, which reflect the deformation of the lumen wall surface under physiological conditions.

[0104] Strain field information is extracted from the real-time deformation data, and a mapping relationship between the strain field information and the nodes of the mesh elements is established. Specifically, the lumen wall surface is divided into 1000 triangular mesh elements, each mesh element has 3 nodes. For the displacement data of each discrete point, the radial basis function interpolation method is used to map it to the corresponding mesh node. For example, for the discrete point located at (10.5mm, 15.2mm, 8.7mm), the displacement value is (0.12mm, -0.08mm, 0.05mm), and the displacement information is transmitted to the surrounding mesh nodes by the interpolation algorithm, so that the mesh deformation is consistent with the actual lumen wall deformation.

[0105] The strain rate tensor of a mesh element is calculated according to the strain field information. For each triangular element, a 3x3 strain rate tensor is constructed based on the displacement differences of its three nodes. This tensor describes the deformation rate of the element in each direction. The strain rate tensor is decomposed into eigenvalues λ1, λ2, λ3 and corresponding eigenvectors v1, v2, v3. The eigenvalues represent the principal strain rates, and the eigenvectors represent the principal strain directions. In a practical case, the principal strain rates of a certain mesh element are 0.15, 0.08, and -0.05, respectively, and the corresponding principal strain directions are (0.866, 0.5, 0), (-0.5, 0.866, 0), and (0, 0, 1), respectively.

[0106] The deformation control parameters of a mesh element are constructed based on the principal strain directions and the principal strain rates. The deformation control parameters include stretch coefficients and rotation angles. The stretch coefficients are proportional to the principal strain rates and are set to 1.2 times the principal strain rates, i.e., 0.18, 0.096, and -0.06 in this example. The rotation angles are determined by aligning the coordinate system of the initial mesh with the principal strain directions. For the above example, the rotation angles are 30°, 120°, and 0°, respectively.

[0107] The node coordinates of a mesh element are adjusted according to the deformation control parameters. For each mesh node, the displacement increment in each principal strain direction is calculated. The displacement increment is equal to the distance of the node to the center of the element multiplied by the stretch coefficient of the corresponding direction. For example, for a node 2 mm away from the center of the element, the displacement increment in the first principal strain direction is 2 mm x 0.18 = 0.36 mm. The node coordinates are updated by considering the displacement increments in the three principal strain directions, and the deformation of the mesh element is realized.

[0108] The strain energy gradient between adjacent mesh elements is calculated after deformation. The strain energy is proportional to the square of the strain rate, and for each element, the strain energy can be represented as the weighted average of the sum of the squares of the three principal strain rates. The strain energy difference between two adjacent elements divided by their center distance is the strain energy gradient. The gradient threshold is set to 0.05 J / mm 3 If the calculated strain energy gradient exceeds this threshold, the node position iterative optimization is triggered.

[0109] The node position iterative optimization employs a gradient descent method to reduce the strain energy difference between adjacent elements by fine-tuning the grid node positions. In each iteration, the nodes are moved towards the element with lower strain energy, and the moving distance is proportional to the strain energy gradient, with a proportionality coefficient of 0.1. For example, if the strain energies of two adjacent elements are 0.2 J and 0.3 J, respectively, and the center distance is 2 mm, the strain energy gradient is 0.05 J / mm, which is exactly equal to the threshold value. If it exceeds the threshold value, such as 0.06 J / mm, the shared node needs to be moved 0.006 mm towards the element with lower strain energy, and the strain energy gradient is recalculated until it is less than or equal to the threshold value of 0.05 J / mm 3 .

[0110] After optimization, a stress interpolation function is established to calculate the stress field distribution on the lumen wall surface. The order of the stress interpolation function is related to the strain rate of the grid element, and the specific relationship is as follows: when the maximum absolute value of the principal strain rate is less than 0.1, a first-order interpolation function is used; when it is between 0.1 and 0.2, a second-order interpolation function is used; and when it is greater than 0.2, a third-order interpolation function is used. In the above example, the maximum absolute value of the principal strain rate is 0.15, so a second-order interpolation function is used.

[0111] The coefficients of the stress interpolation function are determined by satisfying the balance conditions at the grid nodes. For example, for a grid element using a second-order interpolation function, its stress distribution can be represented as a quadratic function of position, containing 6 unknown coefficients. These coefficients are determined by solving a linear equation system composed of node balance equations. The calculation results show that the maximum principal stress on the convex part of the lumen wall surface can reach 0.75 MPa, while in the flat area, the stress value is about 0.3 MPa.

[0112] As Figure 3The distribution of the error of the stress field calculation of the luminal wall under different node optimization strategies is shown. The horizontal axis represents the position along the luminal wall (mm), and the vertical axis represents the calculation error (%). The present technical solution (diamond mark) shows the error variation trend from the initial position of the wall (0 mm) to the distal position (30 mm), with an error range of 2.3%~15.1% and an average error of 8.2%. In comparison, the traditional grid optimization method (square mark) performs similarly in the initial area of the wall, but the error increases sharply in the middle and distal areas (> 15 mm), with a maximum error of 30.5% and an average error of 21.5%. The grid optimization method based on gradient (circle mark) performs between the two, with an error range of 4.7%~22.3% and an average error of 14.8%. The present technical solution performs significantly better in the strain rate sensitive area (15~20 mm), with an error controlled within 10%, which is 20.3 percentage points lower than the traditional method. In particular, in the distal area of the wall (25~30 mm), the error of the present technical solution is stable at about 2.3%, which is reduced by 13.2% compared with the traditional method. This advantage is due to the fact that the present technical solution adopts a strategy that the strain rate is positively correlated with the order of the interpolation function, and the optimal distribution of the node position is achieved through strain energy gradient optimization. Compared with commonly used adaptive grid optimization techniques (such as Spring-based smoothing) and grid refinement methods based on error estimation, the present technical solution pays more attention to the internal relationship between the strain field information and the grid deformation, and achieves more accurate stress field distribution calculation.

[0113] In the present embodiment, by extracting deformation data online and constructing a mapping of the strain field and the grid nodes, the structural deformation can be reflected in real time, and high-speed dynamic monitoring can be achieved. Based on the eigenvalue decomposition of the strain rate tensor, the principal strain direction and the principal strain rate are obtained, and the node position is adjusted after introducing the deformation control parameter, which significantly reduces the grid distortion and interpolation error. By calculating the strain energy gradient between adjacent elements, the iteration optimization is automatically triggered when the threshold is exceeded, so that the grid distribution is adapted to the stress concentration area, and the calculation stability and convergence speed are ensured. The order of the stress interpolation function constructed can be automatically adjusted according to the strain rate, higher order interpolation is used in the high strain area, and the order is simplified in the low strain area, so that the calculation efficiency and result accuracy are considered. The optimized node distribution and adaptive interpolation can finely depict the stress gradient change of the luminal wall, and provide high reliability data for structural strength evaluation and life prediction.

[0114] In an alternative embodiment, constructing a deformation control parameter of a grid element based on the principal strain direction and the principal strain rate comprises:

[0115] The deformation control parameters include a stretching coefficient along the principal strain direction and a compression coefficient perpendicular to the principal strain direction; a ratio of a principal strain rate of the grid cell to a preset reference strain rate is calculated, and the ratio is taken as a stretching coefficient reference value; the stretching coefficient along the principal strain direction is determined according to the stretching coefficient reference value, and the stretching coefficient non-linearly increases with an increase of the ratio;

[0116] A deviation of an aspect ratio of the grid cell from a preset shape factor is calculated, and the compression coefficient perpendicular to the principal strain direction is determined based on the deviation, and the compression coefficient is used to constrain grid distortion;

[0117] The node coordinates of the grid cell are adjusted according to the deformation control parameters, including:

[0118] The node coordinates of the grid cell are projected to a local coordinate system with the principal strain direction as a reference, in which the displacement amount of the node along the principal strain direction is proportional to the stretching coefficient, and the displacement amount of the node perpendicular to the principal strain direction is proportional to the compression coefficient;

[0119] The node coordinates of the grid cell in the global coordinate system are updated based on the displacement amount.

[0120] In a specific embodiment, in the process of real-time positioning and puncture trajectory control in a breast ductoscopy operation, deformation control parameter calculation is a key step to ensure the stability and accuracy of grid deformation. The deformation control parameters mainly include a stretching coefficient along the principal strain direction and a compression coefficient perpendicular to the principal strain direction. For the grid model of the breast tissue, a tetrahedral grid is constructed according to the real-time collected point cloud data, and the principal strain direction and the strain rate of each grid cell need to be calculated. The principal strain direction is obtained by analyzing the deformation gradient tensor of the grid cell, specifically the direction of the eigenvector corresponding to the maximum eigenvalue of the deformation gradient tensor. The principal strain rate of the grid cell is calculated as the length ratio of the grid cell in the principal strain direction before and after deformation, which reflects the stretching or compression degree of the tissue along the direction. The preset reference strain rate is the strain rate of the breast tissue in a normal physiological state, and is usually taken as a value between 1.05 and 1.15, which is obtained based on a large amount of clinical data statistics.

[0121] The ratio of the principal strain rate of the mesh element to the preset reference strain rate is calculated as the reference value of the stretch coefficient. For example, if the principal strain rate of a mesh element is 0.18 and the preset reference strain rate is 0.12, the ratio is 1.5. The stretch coefficient has a nonlinear relationship with the ratio, and the final stretch coefficient is determined by using an S-shaped function mapping relationship. When the ratio is less than 0.8, the stretch coefficient tends to 0.9; when the ratio is between 0.8 and 1.2, the stretch coefficient is approximately linearly related to the ratio; when the ratio is greater than 1.2, the stretch coefficient gradually slows down, and when the ratio reaches 2.0, the stretch coefficient tends to the upper limit value 1.6. This nonlinear mapping relationship can effectively prevent mesh distortion caused by excessive stretching. In practical applications, for the fibrous tissue region around the milk duct, the ratio is usually between 1.3 and 1.5, and the corresponding stretch coefficient is about 1.35 to 1.45; and for the adipose tissue region, the ratio is usually between 0.9 and 1.1, and the corresponding stretch coefficient is about 0.95 to 1.05.

[0122] The aspect ratio of the mesh element is an important indicator for evaluating the quality of the mesh, which is defined as the ratio of the length of the mesh element in the principal strain direction to the width perpendicular to the principal strain direction. The preset shape factor is the ideal aspect ratio of the mesh element, and for the tetrahedral mesh of the tissue around the milk duct, the value is usually set to be between 1.2 and 1.5. The deviation of the actual aspect ratio of the mesh element from the preset shape factor is calculated, and the greater the deviation, the more serious the mesh distortion. The compression coefficient is determined based on the deviation and is used to constrain the deformation of the mesh perpendicular to the principal strain direction. When the deviation is negative, i.e., the aspect ratio is less than the preset shape factor, the compression coefficient linearly increases, and the maximum value is 1.2; when the deviation is positive and does not exceed 0.3, the compression coefficient linearly decreases; when the deviation exceeds 0.3, the compression coefficient rapidly decreases, and the minimum value is 0.7, to prevent excessive stretching of the mesh. For example, the aspect ratio of a mesh element is 1.8 and the preset shape factor is 1.5, so the deviation is 0.3, and the corresponding compression coefficient is about 0.85.

[0123] When adjusting the node coordinates of the mesh element according to the calculated deformation control parameters, a local coordinate system with the principal strain direction as the reference is needed to be established. The x-axis of the local coordinate system is along the principal strain direction, and the y-axis is perpendicular to the principal strain direction. The node coordinates of the mesh element are converted from the global coordinate system to the local coordinate system, and the conversion matrix is composed of the principal strain direction vector and its orthogonal vector. In the local coordinate system, the displacement amount of the node along the x-axis (the principal strain direction) is proportional to the stretching coefficient, and the calculation method is to multiply the original x-coordinate of the node by the stretching coefficient and then subtract the original x-coordinate. For example, the x-coordinate of a node in the local coordinate system is 2.5 mm, and the stretching coefficient is 1.4, then the x-direction displacement amount is 2.5 x 1.4 - 2.5 = 1.0 mm. The displacement amount of the node along the y-axis (perpendicular to the principal strain direction) is proportional to the compression coefficient, and the calculation method is to multiply the original y-coordinate of the node by the compression coefficient and then subtract the original y-coordinate. For example, the y-coordinate of a node in the local coordinate system is -1.2 mm, and the compression coefficient is 0.85, then the y-direction displacement amount is -1.2 x 0.85 + 1.2 = -0.18 mm.

[0124] The node coordinates of the mesh element in the local coordinate system are updated based on the calculated displacement amount, and then the updated node coordinates are converted back to the global coordinate system through the inverse matrix of the coordinate conversion matrix. In the breast ductoscopy surgery, the number of typical mesh elements of the breast duct surrounding tissue is 5000 to 10000, and for each mesh element, the above deformation control parameter calculation and node coordinate updating operations need to be performed. Through GPU accelerated calculation, the entire mesh deformation process can be completed within 50 milliseconds, meeting the real-time control requirements of the surgery.

[0125] In the prior art, linear elastic model or particle-spring system is mainly used to process breast tissue deformation. These methods have obvious limitations: linear elastic model cannot accurately describe the nonlinear deformation behavior under large deformation; particle-spring system is difficult to maintain volume constraint, and is easy to cause unrealistic volume loss; in addition, the existing methods mostly use uniform deformation parameters, ignoring the deformation characteristic differences of different tissue types and regions, and cannot accurately simulate the heterogeneous deformation of breast tissue. In terms of mesh quality control, traditional methods mainly rely on simple stiffness coefficient or Laplace smoothing, lack of accurate control of local shape of mesh, and are easy to cause excessive smoothing or mesh degradation.

[0126] The method of the embodiment proposes a deformation control parameter calculation method based on principal strain analysis. Compared with the prior art, the concepts of principal strain direction and principal strain rate are introduced, so that the deformation parameters are directly related to the actual deformation characteristics of the tissue; a nonlinear mapping relationship is used to determine the stretching coefficient, effectively avoiding the distortion problem of the linear model in the case of large deformation; a compression coefficient is designed based on the grid shape factor, realizing accurate control of the grid quality; anisotropic deformation is realized through local coordinate system transformation, accurately simulating the deformation difference of breast tissue in different directions. While ensuring the real-time performance of deformation calculation, the accuracy and stability of soft tissue deformation simulation are improved. By introducing the deformation control parameters and the local coordinate transformation mechanism, the method of the embodiment can realize accurate control of the deformation of the breast tissue in a complex surgical environment.

[0127] In an optional embodiment, according to the stress field distribution, a puncture safety coefficient is calculated through multi-layer sampling points, an optimal puncture trajectory is determined based on the puncture safety coefficient, and the movement tracking control of the puncture instrument is realized in combination with a speed control point sequence, including:

[0128] According to the stress field distribution, multi-layer sampling points are constructed in the surgical space, and each layer of sampling points is connected to form a sampling plane; a connection relationship is established between adjacent sampling planes, and a puncture safety coefficient between sampling points is calculated according to the connection relationship, wherein the puncture safety coefficient is inversely proportional to the stress field value at the sampling point;

[0129] A starting point and a target point of puncture are set, a virtual potential field distribution map based on the puncture safety coefficient is constructed, the direction of the resultant force of the attractive force and the repulsive force is calculated in the virtual potential field distribution map, and the optimal puncture path is calculated through adaptive step gradient iteration;

[0130] Speed control points are set along the optimal puncture path, and the movement speed of the puncture instrument at each speed control point is determined according to the change of the lumen curvature;

[0131] A tangent vector and a normal vector are obtained at each speed control point, and the movement posture parameters of the puncture instrument are calculated;

[0132] Real-time movement parameters of the puncture instrument are collected, the real-time movement parameters are compared with the movement posture parameters, and a deviation compensation signal is generated;

[0133] According to the deviation compensation signal, the driving output of the puncture instrument is adjusted to realize continuous movement control of the puncture instrument along the optimal puncture path.

[0134] In a specific embodiment, when constructing the multi-layer sampling points in the surgical space, a plurality of planes perpendicular to the puncture direction can be set at intervals of 5 mm in the puncture space, 100x100 sampling points are uniformly distributed on each plane, and a total of 20 planes are set. The sampling points between adjacent planes establish a connection relationship according to the eight-neighborhood principle, that is, each sampling point establishes a connection with the eight nearest points of the adjacent plane. The puncture safety coefficient C between the sampling point (i, j, k) and the sampling point (i+1, j+1, k+1) can be calculated by the inverse ratio of the stress field value. For example, when the stress value at the sampling point (i, j, k) is 5 MPa and the stress value at the sampling point (i+1, j+1, k+1) is 8 MPa, the safety coefficient of the path between the two points can be calculated as C = 1 / (5+8) = 0.0769. The higher the safety coefficient, the safer the path.

[0135] When setting the starting point and target point of the puncture, it can be determined according to the doctor's mark. For example, the starting point is set as the skin surface coordinate (50, 50, 0), and the target point is set as the lesion center point coordinate (80, 75, 85). When constructing the virtual potential field distribution map, the safety coefficient of each sampling point and the distance from the target point are considered comprehensively. The target point has the lowest potential energy value 0, and the high-risk area (such as the area near the large blood vessel, the area with a stress value greater than 10 MPa) sets a higher potential energy value, such as 100. For general tissue areas, the potential energy value can be linearly mapped to the range of 10-50 according to the stress field value.

[0136] When calculating the direction of the resultant force in the virtual potential field, the attractive force direction always points to the target point, the attractive force size is proportional to the distance from the current point to the target point, and the attractive force coefficient is set to 0.5. The repulsive force direction is perpendicular to the high stress gradient direction, and the repulsive force size is proportional to the stress field gradient value, and the repulsive force coefficient is set to 2.0. In each iteration step, the next moving direction is determined according to the direction of the resultant force of the attractive force and the repulsive force.

[0137] When calculating the adaptive step gradient iteration, the initial step is set to 2 mm, and the step is dynamically adjusted according to the local stress field gradient during the iteration process. When approaching the high-risk area (stress gradient greater than 1.5 MPa / mm), the step is reduced to 0.5 mm; in the safe area (stress gradient less than 0.3 MPa / mm), the step is increased to 3 mm. Iterative calculation continues until the starting point reaches the target point, or the maximum iteration number 500 is reached. Through this calculation, an optimal puncture path connecting the starting point and the target point while avoiding the high-risk area can be obtained, which is composed of a series of space point coordinates.

[0138] After the optimal puncture path is obtained, speed control points are evenly set along the path. When the total length of the path is 105 mm, one control point can be set every 5 mm, and a total of 21 control points are obtained. The motion speed of each control point is determined according to the change of the curvature of the lumen, and the curvature is calculated based on the curvature radius formed by the adjacent three points. For example, for the i-th point in the control point sequence, its curvature K can be approximately calculated by the coordinate difference of the adjacent points before and after it. When the curvature at the control point exceeds 0.2 / mm (corresponding to a curvature radius of less than 5 mm), the puncture speed is reduced to 30% of the basic speed, about 1.5 mm / s; when the curvature is between 0.1-0.2 / mm, the speed is set to 60% of the basic speed, about 3 mm / s; when the curvature is less than 0.1 / mm, the basic speed of 5 mm / s is maintained.

[0139] At each speed control point, a tangent vector and a normal vector are calculated based on the adjacent control points. The tangent vector points to the next control point, and the normal vector is determined by the normal of the plane formed by the previous point and the next point. Taking the control point (65, 60, 35) as an example, its previous and next control points are (63, 58, 30) and (68, 62, 40) respectively. According to the three points, the tangent vector at this point is (5, 4, 10), and after unitization, it is (0.41, 0.33, 0.82), and the normal vector is (-0.82, 0.57, 0.03). The posture parameters of the puncture instrument include position coordinates, direction angles (pitch angle 23°, yaw angle 38°) and speed 3 mm / s.

[0140] In the actual puncture process, the real-time motion parameters of the puncture instrument are collected by the optical tracking system every 100 ms, including spatial position, attitude angle and motion speed. These real-time parameters are compared with the predetermined motion posture parameters to generate a deviation compensation signal. For example, when the actual position of the puncture instrument is (65.3, 59.7, 35.2), which deviates from the ideal position (65, 60, 35), the position deviation vector (0.3, -0.3, 0.2) is calculated; when the actual attitude angle is (pitch angle 24.5°, yaw angle 37.2°), the angle deviation is (1.5°, -0.8°).

[0141] The driving output of the puncture instrument is adjusted according to the deviation compensation signal, and a PID control algorithm is used to realize accurate compensation. The PID parameters of position control are set to P=2.5, I=0.3, and D=0.8; the PID parameters of attitude angle control are set to P=1.8, I=0.2, and D=0.5. When the position deviation vector (0.3, -0.3, 0.2) is detected, the compensation driving force is calculated as (-0.75, 0.75, -0.5) N; when the angle deviation (1.5°, -0.8°) is detected, the compensation torque is calculated as (-2.7, 1.44) N·mm. Through this real-time closed-loop control, the stable movement of the puncture instrument along the optimal path is ensured, and the positioning accuracy can reach ±0.5 mm, and the angle accuracy can reach ±1°.

[0142] In an alternative embodiment, the starting point and target point of puncture are set, a virtual potential field distribution map based on the puncture safety coefficient is constructed, the direction of the resultant force of the attractive force and the repulsive force is calculated in the virtual potential field distribution map, and the optimal puncture path is calculated by adaptive step gradient iteration, including:

[0143] The spatial coordinates of the puncture starting point and target point are determined on the sampling plane, a basic potential field function is constructed with the reciprocal of the puncture safety coefficient, a local potential field distribution is established at each sampling point, and the local potential fields of adjacent sampling planes are interpolated and fitted to form a continuous virtual potential field distribution map;

[0144] An attractive force field is established with the target point as the attractive force source, the strength of the attractive force field decays with increasing distance, a repulsive force field is established with the sampling points whose puncture safety coefficients are less than a preset safety threshold as the repulsive force source, the strength of the repulsive force field is inversely proportional to the puncture safety coefficient, the vector sum of the attractive force and the repulsive force is superimposed at the calculation points of the virtual potential field distribution map, and the direction of the resultant force is determined;

[0145] An adaptive step gradient descent algorithm is used, the starting point is taken as the initial search position, iterative calculation is performed along the direction of the resultant force until the position error between the search position and the target point is less than a preset error threshold, the optimal puncture path is obtained, and the optimal puncture path is smoothed to ensure continuous curvature.

[0146] In a specific embodiment, in the process of real-time positioning and trajectory control in the ductoscopy surgery, the analysis of the sampling plane is the basis for determining the optimal puncture path. The sampling plane refers to a plane perpendicular to the center line of the milk duct set every 1 mm along the direction of the milk duct. There are 64 sampling points uniformly distributed on each sampling plane for evaluating the puncture safety at different positions on the plane. The puncture starting point is usually selected at the position closest to the target lesion on the surface of the nipple, and its spatial coordinates are determined by preoperative three-dimensional ultrasound or MRI images. The target point is the center position of the lesion that needs to be reached, which is specified by the doctor when planning preoperatively. For example, in a typical ductoscopy surgery, the puncture starting point coordinates can be (5.2mm, 8.7mm, 0.0mm), and the target point coordinates can be (42.6mm, 35.8mm, 28.3mm), which are local coordinate systems established relative to the nipple tip.

[0147] The puncture safety coefficient is a key indicator representing the degree of puncture safety, with a value range of 0 to 1, and the larger the value, the higher the puncture safety. The calculation of the puncture safety coefficient considers the distance between the sampling point and important anatomical structures (such as blood vessels, nerves, etc.), the density of the surrounding tissue, and the elastic properties of the tissue at that point. The basic potential field function is constructed based on the reciprocal of the puncture safety coefficient, i.e., the safety coefficient s, and the basic potential field value is 1 / s. The purpose of this design is to make the area with low safety coefficient have high potential energy, forming a "potential energy mountain", and the area with high safety coefficient has low potential energy, forming a "potential energy valley". For example, the safety coefficient of a point is 0.2, and its basic potential field value is 5; the safety coefficient of a point is 0.8, and its basic potential field value is 1.25. A local potential field distribution is established at each sampling point, and the local potential field distribution is usually described by a radial basis function with an influence range of 2mm. For the area between adjacent sampling planes, a cubic spline interpolation method is used to fit the local potential field, forming a continuous virtual potential field distribution map, which ensures the continuity and smoothness of the potential field in the entire three-dimensional space.

[0148] The design of the attractive field and the repulsive field is the core of the virtual potential field method. The target point is the attractive source, and the attractive field is established to make the puncture path develop towards the target point. The attractive field strength decays with the increase of the distance from the target point, and a quadratic decay model is adopted, i.e., the attractive strength at a point with a distance of d is k / (d 2+ ε), where k is an attractive force coefficient, usually set to 10, and ε is a small constant to prevent the denominator from being zero, taking a value of 0.01. When the distance between the sampling point and the target point is less than 5 mm, the attractive force strength remains at a maximum value; when the distance is greater than 20 mm, the attractive force strength approaches an asymptote. The repulsive force field is generated by sampling points with a puncture safety coefficient less than a preset safety threshold, which is usually set to 0.5. The strength of the repulsive force field is inversely proportional to the puncture safety coefficient and is calculated as m x (1 - s) / s, where m is a repulsive force coefficient, usually set to 5, and s is the puncture safety coefficient. For example, a point with a safety coefficient of 0.3 has a repulsive force strength of 11.67, and a point with a safety coefficient of 0.4 has a repulsive force strength of 7.5. The range of the repulsive force is usually limited to within 3 mm to avoid excessive influence of distant dangerous areas on path planning. At each calculation point in the virtual potential field distribution map, all the attractive and repulsive forces acting on the point are vector superimposed to obtain the direction of the resultant force, which indicates the locally optimal puncture advancement direction.

[0149] The gradient descent algorithm with adaptive step size is used to find the optimal puncture path in a complex virtual potential field. The puncture starting point is taken as the initial search position, and the search position is iteratively updated along the direction of the resultant force. The adaptive step size means that the step size is not a fixed value, but is dynamically adjusted according to the size of the potential field gradient at the current position, and is calculated as α / (1 + β x |G|), where α is the basic step size, usually set to 0.5 mm, β is the adjustment coefficient, taking a value of 2, and |G| is the modulus of the potential field gradient. In areas with large potential field gradients, the step size is automatically reduced to improve search accuracy; in flat areas with small potential field gradients, the step size is automatically increased to speed up convergence. For example, when the modulus of the potential field gradient is 0.2, the step size is about 0.45 mm; when the modulus of the potential field gradient is 2.0, the step size is about 0.2 mm. Iterative calculation continues until the Euclidean distance between the current search position and the target point is less than a preset error threshold, which is usually set to 0.8 mm. In practical applications, the number of iterations is generally between 100 and 200, and the entire path planning process can be completed within 200 milliseconds, meeting the real-time operation requirements.

[0150] After obtaining the optimal puncture path, smoothing processing is required to ensure curvature continuity and avoid sharp turns in the path. The smoothing processing uses a cubic B-spline curve fitting method, taking the original path points as control points to generate a smooth curve with continuous curvature. During smoothing, the maximum deviation of the generated curve from the original path should not exceed 0.5 mm, and the smoothed path should still meet the safety requirements. After smoothing, the final path is discretized into an equidistant point sequence, with a point spacing of 0.5 mm, to facilitate subsequent puncture mechanism execution. For a typical puncture path of 40 mm in length, the number of path points after smoothing is about 80, and the three-dimensional coordinates of these points constitute the complete puncture trajectory data.

[0151] The method of the embodiment is mainly based on the artificial potential field method in the field of robot path planning and the surgical navigation technology under the guidance of medical images. The existing technology mainly adopts the geometric constraint method based on anatomical structure segmentation or the A algorithm based on voxels for puncture path planning in the ductoscopy surgery. The geometric constraint method sets dangerous structures such as blood vessels and nerves as forbidden zones by pre-segmenting these structures, and plans the shortest path around these forbidden zones, but this method is complex in calculation and difficult to update the path in real time. The A algorithm based on voxels discretizes the patient's anatomical space into a voxel grid, and searches for the optimal path on the grid, but there is a trade-off between voxel size and search accuracy and computational efficiency, and it is difficult to adapt to tissue deformation. The existing methods are mostly based on static preoperative images, and cannot adapt to the dynamic deformation of tissues during the surgery, increasing the risk of puncture deviation.

[0152] The method of the embodiment proposes an adaptive puncture path planning method based on virtual potential field. The puncture safety coefficient is introduced into the potential field construction process, so that the path planning is directly related to the safety of the surgery; the continuous potential field distribution is used instead of the discrete forbidden zone representation, which improves the smoothness and naturalness of the path; the gradient descent algorithm with adaptive step size is introduced to balance the computational efficiency and path accuracy; the potential field distribution is dynamically updated combined with the real-time acquired tissue deformation information, so that the path planning can adapt to the tissue deformation. Under the condition of tissue deformation, a safe and smooth puncture path is planned in real time, which avoids important anatomical structures and accurately reaches the target lesion. By introducing the virtual potential field based on the safety coefficient and the adaptive search strategy, the method of the embodiment can significantly improve the safety and accuracy of the puncture path while ensuring the computational efficiency.

[0153] As Figure 4As shown, a virtual potential field-based puncture path planning effect diagram is displayed, which intuitively presents the advantages of the technical scheme compared with the prior art. In the figure, the virtual potential field distribution is shown in the form of a heat map, and the color from green to red represents the potential field strength from low to high (the safety factor from high to low). The dangerous area is marked with a square rotated by 45 degrees, and the safety factor value is less than 0.4. The starting point is located at coordinates (1, 1), and the target point is located at coordinates (18, 18). The path (solid square mark) planned by the technical scheme successfully avoids three main dangerous areas: high potential energy areas near (5, 10), (12, 6) and (15, 15). From the path diagram, it can be observed that the path generated by the technical scheme is smoother, and the curvature is continuous. The turning radius at position (7, 7) is about 2.5 mm, while the turning of the RRT algorithm (triangular mark dotted line) at the same position presents a sharp corner, and the turning radius is about 0.8 mm. In addition, the total length of the path of the technical scheme is 24.3 mm, while the path length of the RRT algorithm is 26.7 mm, and the path length of the A algorithm (circular mark dotted line) is 29.5 mm, which significantly improves the path length. In terms of safety, the minimum safety factor of the path of the technical scheme is 0.52 (located at coordinates (8, 9)), which is significantly higher than 0.38 of the RRT algorithm and 0.35 of the A algorithm, which reflects the technical advantage of the technical scheme in realizing path safety and smoothness through the superposition of gravitational field and repulsive field.

[0154] The real-time positioning and puncture trajectory control system in the mammary ductoscopy surgery of the embodiment of the application comprises:

[0155] A first unit is configured to collect an endoscopic image sequence of a mammary ductoscopy probe, acquire spatial positioning data through an electromagnetic position sensor on the probe, time sequence register the endoscopic image sequence and the spatial positioning data, and establish a real-time coordinate system of a surgery space;

[0156] A second unit is configured to process the endoscopic image sequence through a light field reconstruction algorithm based on the real-time coordinate system, acquire three-dimensional depth information of a lumen inner wall, determine real-time deformation parameters of a lumen channel in combination with the spatial positioning data of the electromagnetic position sensor, and determine real-time deformation parameters of a lumen channel in combination with the spatial positioning data of the electromagnetic position sensor.

[0157] A third unit is configured to perform adaptive deformation on the lumen wall surface by using a grid form adjustment method, control node optimization by using a strain energy gradient, and calculate stress field distribution of the lumen wall surface according to the real-time deformation data.

[0158] A fourth unit is configured to calculate a puncture safety factor through a plurality of sampling points according to the stress field distribution, determine an optimal puncture trajectory based on the puncture safety factor, and realize motion tracking control of a puncture instrument in combination with a speed control point sequence.

[0159] A fifth unit is configured to acquire deformation signals of the puncture instrument in real time by using a fiber strain sensing array, and calculate real-time force data of the puncture instrument.

[0160] A sixth unit is configured to analyze deviation between the actual motion track of the puncture instrument and the puncture track according to the real-time force data, to adjust the compensation parameter in real time, and to realize precise track tracking control of the puncture instrument.

[0161] In a third aspect, the present application provides an electronic device, comprising:

[0162] a processor;

[0163] a memory for storing processor-executable instructions;

[0164] The processor is configured to call the instructions stored in the memory to execute the method described above.

[0165] In a fourth aspect, the present application provides a computer-readable storage medium, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method described above.

[0166] The present application can be a method, device, system and / or computer program product. The computer program product can include a computer-readable storage medium having stored thereon computer-readable program instructions that, when executed by a computer, cause the computer to carry out various aspects of the present application.

[0167] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A real-time positioning and puncture trajectory control system for ductoscopy, characterized in that, include: The first unit is used to acquire endoscopic image sequences from the ductoscope probe and obtain spatial positioning data through the electromagnetic position sensor on the probe. The endoscopic image sequences and spatial positioning data are then time-series registered to establish a real-time coordinate system for the surgical space. The second unit is used to set up an optical field sensor array in the real-time coordinate system, collect multi-view light information of the inner wall of the tube through the optical field sensor array, and perform optical field reconstruction processing on the multi-view light information to obtain the main ray coordinates and micro ray coordinates. The spatial reprojection error of the ray is calculated based on the coordinates of the principal ray and the coordinates of the micro ray. An optimization objective is constructed by combining the depth smoothing constraint term. The three-dimensional depth information of the inner wall of the tube is obtained by iterative optimization using the gradient descent method. The spatial positioning data of the electromagnetic position sensor and the three-dimensional depth information are transformed by coordinate transformation to generate three-dimensional point cloud data characterizing the inner wall of the tube. Feature matching and registration are performed on continuously acquired 3D point cloud data of multiple frames to extract point cloud deformation vectors. Real-time deformation parameters of the lumen channel are obtained by introducing elastic constraints and iterative optimization. The third unit is used to adaptively deform the tube wall using a mesh shape adjustment method, optimize the nodes by controlling the strain energy gradient, and calculate the stress field distribution of the tube wall based on the real-time deformation parameters. The fourth unit is used to calculate the puncture safety factor through multiple sampling points based on the stress field distribution, determine the optimal puncture trajectory based on the puncture safety factor, and realize the motion tracking control of the puncture instrument by combining the speed control point sequence. The fifth unit is used to acquire the deformation signal of the puncture instrument in real time using a fiber optic strain sensor array and calculate the real-time force data of the puncture instrument. The sixth unit is used to analyze the deviation between the actual movement trajectory of the puncture instrument and the puncture trajectory based on the real-time force data, and adjust the compensation parameters in real time to achieve precise trajectory tracking control of the puncture instrument.

2. The system according to claim 1, characterized in that, Feature matching and registration are performed on continuously acquired multi-frame 3D point cloud data to extract point cloud deformation vectors. Real-time deformation parameters of the lumen channel are obtained through elastic constraints and iterative optimization, including: The curvature, normal vector and geometric moment distribution of point cloud data are calculated from multiple frames of continuously acquired 3D point cloud data to obtain feature descriptors that characterize the local geometric properties of each point cloud. Based on the feature descriptor, the K-nearest neighbor search algorithm is used to perform feature matching on the point clouds of adjacent frames to obtain initial matching point pairs. The transformation relationship between the initial matching point pairs is calculated iteratively, and abnormal matching point pairs are removed to obtain the point cloud registration result. Based on the point cloud registration results, the deformation direction and deformation amplitude of the local point cloud are extracted, a point cloud deformation vector is constructed, and gradient analysis is performed on the point cloud deformation vector to obtain the deformation distribution characteristics. Regional deformation constraints are constructed based on the point cloud deformation vector transfer of nearest neighbor pairs, and the corresponding elastic correction values ​​are calculated. The deformation gradient and elastic deformation of adjacent regions of the point cloud are extracted based on the elastic correction values. A deformation parameter mapping model is established using the moving least squares method. The point cloud deformation vector is iteratively optimized in multiple rounds based on the deformation parameter mapping model. The real-time deformation parameters of the cavity channel are determined by calculating the difference value of the point cloud deformation vector in adjacent iteration rounds.

3. The system according to claim 1, characterized in that, An adaptive deformation method for the tube wall is employed using a mesh shape adjustment approach. Strain energy gradient control is used for node optimization. The stress field distribution on the tube wall is calculated based on the real-time deformation parameters, including: Extract strain field information from real-time deformation parameters and establish a mapping relationship between the strain field information and mesh element nodes; The strain rate tensor of the mesh element is calculated based on the strain field information, and the principal strain direction and principal strain rate of the mesh element are obtained by eigenvalue decomposition of the strain rate tensor. Deformation control parameters for the mesh elements are constructed based on the principal strain direction and the principal strain rate; the node coordinates of the mesh elements are adjusted according to the deformation control parameters. Calculate the strain energy gradient between adjacent grid cells after adjustment. When the strain energy gradient exceeds a preset gradient threshold, trigger iterative optimization of node positions until the strain energy gradient is less than or equal to the gradient threshold. A stress interpolation function is established based on the optimized grid node distribution. The order of the stress interpolation function is positively correlated with the strain rate of the corresponding grid element, and the stress field distribution on the tube wall is calculated.

4. The system according to claim 3, characterized in that, The deformation control parameters for constructing mesh elements based on the principal strain direction and the principal strain rate include: The deformation control parameters include the tensile coefficient along the principal strain direction and the compressive coefficient perpendicular to the principal strain direction; the ratio of the principal strain rate of the mesh element to the preset reference strain rate is calculated, and the ratio is used as the reference value of the tensile coefficient. The tensile coefficient along the principal strain direction is determined based on the reference value of the tensile coefficient, and the tensile coefficient increases nonlinearly with the increase of the ratio. The deviation between the aspect ratio of the mesh cells and the preset shape factor is calculated, and the compression coefficient perpendicular to the principal strain direction is determined based on the deviation. The compression coefficient is used to constrain mesh distortion. Adjusting the node coordinates of the mesh elements according to the deformation control parameters includes: The node coordinates of the mesh element are projected onto a local coordinate system based on the principal strain direction. In the local coordinate system, the displacement of the node along the principal strain direction is proportional to the tensile coefficient, and the displacement of the node perpendicular to the principal strain direction is proportional to the compressibility coefficient. The node coordinates of the mesh element in the global coordinate system are updated based on the displacement.

5. The system according to claim 1, characterized in that, Based on the stress field distribution, a puncture safety factor is calculated using multi-layer sampling points. Based on this puncture safety factor, the optimal puncture trajectory is determined. Combined with a speed control point sequence, motion tracking control of the puncture instrument is achieved, including: Based on the stress field distribution, multiple sampling points are constructed in the surgical space, and each layer of sampling points is connected to form a sampling plane; a connection relationship is established between adjacent sampling planes, and a puncture safety factor between sampling points is calculated based on the connection relationship, wherein the puncture safety factor is inversely proportional to the stress field value at the sampling point; Set the starting point and target point of the puncture, construct a virtual potential field distribution map based on the puncture safety factor, calculate the resultant force direction of the attraction and repulsion in the virtual potential field distribution map, and calculate the optimal puncture path through adaptive step size gradient iteration. Speed ​​control points are set along the optimal puncture path, and the movement speed of the puncture instrument at each speed control point is determined according to the change in lumen curvature. At each velocity control point, the tangent vector and normal vector are obtained, and the motion posture parameters of the puncture instrument are calculated. Real-time motion parameters of the puncture instrument are collected, and the real-time motion parameters are compared with the motion posture parameters to generate a deviation compensation signal. The drive output of the puncture instrument is adjusted according to the deviation compensation signal to achieve continuous motion control of the puncture instrument along the optimal puncture path.

6. The system according to claim 5, characterized in that, Set the starting and target points for the puncture, construct a virtual potential field distribution map based on the puncture safety factor, calculate the resultant force direction of attraction and repulsion in the virtual potential field distribution map, and calculate the optimal puncture path through adaptive step-size gradient iteration, including: The spatial coordinates of the puncture start point and target point are determined on the sampling plane. The basic potential field function is constructed using the reciprocal of the puncture safety factor. A local potential field distribution is established at each sampling point. The local potential fields of adjacent sampling planes are interpolated and fitted to form a continuous virtual potential field distribution map. A gravitational field is established with the target point as the gravitational source. The strength of the gravitational field decreases with increasing distance. A repulsive field is established with the sampling point where the puncture safety factor is less than the preset safety threshold as the repulsive source. The strength of the repulsive field is inversely proportional to the puncture safety factor. The vector sum of the gravitational and repulsive forces is superimposed at the calculation point of the virtual potential field distribution map to determine the direction of the resultant force. An adaptive step-size gradient descent algorithm is used to perform iterative calculations along the direction of the resultant force, starting from the initial search position, until the positional error between the search position and the target point is less than a preset error threshold, thereby obtaining the optimal puncture path. The optimal puncture path is then smoothed to ensure curvature continuity.

Images