Interventional ultrasound image real-time analysis method based on feature point extraction and elastic registration

Abstract

The application discloses an interventional ultrasound image real-time analysis method based on feature point extraction and elastic registration, and belongs to the technical field of interventional ultrasound navigation and medical image processing. The method comprises the following steps: acquiring reference body data and generating reference safety boundary body data; continuously receiving image and pose data of an interventional ultrasound probe, and performing preprocessing and space-time calibration; extracting feature points of an anatomical structure and an interventional instrument and constructing feature point graph information; performing rigid registration and elastic registration; performing interventional instrument positioning and lesion ablation state analysis, and generating monitoring indexes and risk prompt information. The application aims at the problem that in the prior art, soft tissue deformation, unclear lesion boundary and unstable instrument positioning during operation lead to difficulty in realizing whole-process real-time quantitative monitoring. The real-time compensation for intraoperative deformation is realized through elastic registration driven by a feature point graph constraint, so that the real-time fusion display and quantitative analysis of the interventional ultrasound are realized without performing intraoperative three-dimensional reconstruction, and the safety is improved.

Description

Technical Field

[0001] This invention belongs to the field of interventional ultrasound navigation and medical image processing technology, specifically involving a real-time analysis method for interventional ultrasound images based on feature point extraction and elastic registration. It is applicable to the real-time analysis of soft tissue deformation compensation, device positioning and efficacy evaluation during interventional procedures such as ablation, puncture, biopsy and catheter placement. Background Technology

[0002] Interventional ultrasound is widely used in clinical settings such as radiofrequency ablation, microwave ablation, puncture biopsy, drainage catheter placement, and nerve blocks. Ultrasound imaging offers advantages such as real-time performance, low cost, and no radiation. However, during interventional procedures, factors such as respiratory movements, changes in body position, water isolation, probe pressure, and vaporization caused by thermal ablation can significantly deform and displace the target organ and lesion. This results in unclear lesion boundaries on ultrasound images, making it difficult to stably locate the lesion and its safe boundaries. Physicians often rely on experience for judgment, leading to significant subjective errors.

[0003] Existing technologies include various schemes related to ultrasound image registration or target tracking. For example, local offsets are estimated through image patch matching, feature point matching, or contour matching to achieve sampling gate or target region tracking. However, these schemes are mostly used for translation compensation within two-dimensional images and are difficult to effectively model three-dimensional soft tissue deformation. Other technologies achieve flexible registration of preoperative body data and intraoperative ultrasound sequences through single registration combined with multiple secondary fine registrations. However, these often rely on strong segmentation quality and lack assessment and degradation strategies for registration confidence, making them prone to registration drift when instruments are obstructed or echo quality deteriorates. Still other technologies achieve dynamic three-dimensional reconstruction through point cloud-level rigid / elastic registration, but this usually requires additional three-dimensional reconstruction or point cloud generation processes, making it difficult to meet the real-time and stability requirements of interventional ultrasound.

[0004] Therefore, there is an urgent need for a real-time interventional ultrasound image analysis method that can drive elastic registration using stable and traceable feature point information in interventional ultrasound sequences without intraoperative three-dimensional reconstruction, thereby achieving real-time compensation for organ deformation and outputting quantifiable monitoring indicators. Summary of the Invention

[0005] This application provides a real-time analysis method for interventional ultrasound images based on feature point extraction and elastic registration. It aims to address the problems in existing technologies where soft tissue deformation during intervention makes it difficult to stably locate lesions and safety boundaries, and where instrument positioning relies on experience and lacks quantitative monitoring standards throughout the entire process. The method achieves rigid / elastic joint registration driven by feature point map constraints, generating deformed reference volume data and safety boundaries without intraoperative 3D reconstruction. It further outputs instrument positioning results, ablation coverage, and risk warning information. The method includes: acquiring reference volume data and generating reference safety boundary volume data; continuously receiving two-dimensional ultrasound image frames output by the interventional ultrasound probe and probe pose data synchronized with the time of the two-dimensional ultrasound image frames; and performing preprocessing and timing analysis on the two-dimensional ultrasound image frames. The system performs empty calibration to generate a spatiotemporally consistent interventional ultrasound image sequence. It continuously reads the interventional ultrasound image sequence, extracts feature points of anatomical structures and interventional instruments, constructs feature point map information, and generates a feature point map information sequence with confidence. It continuously reads the feature point map information sequence and reference body data, performs rigid and elastic registration on the reference body data based on the feature point map information, generates a rigid transformation matrix and a three-dimensional elastic deformation field, and generates deformed reference body data and deformed safety boundaries. It continuously reads real-time two-dimensional ultrasound images and deformed reference body data, performs interventional instrument positioning and lesion ablation status analysis, generates monitoring indicators and risk warning information, and outputs them. It continuously monitors registration confidence and abnormal states, and performs downgrading and reinitialization control when abnormal states are identified.

[0006] One or more technical solutions provided in this application have at least the following technical effects or advantages: This application provides a real-time analysis method for interventional ultrasound images based on feature point extraction and elastic registration, relating to the fields of interventional ultrasound navigation and medical image processing. By extracting feature points of anatomical structures and interventional instruments from the interventional ultrasound sequence and constructing feature point maps, the method introduces feature point map constraints into a rigid / elastic joint registration process. This enables real-time compensation for intraoperative deformation in scenarios involving soft tissue deformation caused by probe compression, respiration, and thermal ablation. Furthermore, it achieves dynamic updating of reference body data and real-time mapping of safety boundaries without intraoperative 3D reconstruction, outputting instrument positioning results, ablation coverage, and risk warning information. This improves the stability, interpretability, and quantification of the entire interventional procedure monitoring process, reduces the risk of misjudgment, and enhances the safety and effectiveness of ablation / puncture procedures. Attached Figure Description

[0007] To more clearly illustrate the technical solutions in the embodiments of this application, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0008] Figure 1 This is a schematic diagram of a real-time analysis method for interventional ultrasound images based on feature point extraction and elastic registration, provided as an embodiment of this application. Detailed Implementation

[0009] This application provides a real-time analysis method for interventional ultrasound images based on feature point extraction and elastic registration, which solves the technical problems in the prior art where soft tissue deformation during intervention makes it difficult to stably locate lesions and safety boundaries, and where instrument positioning relies on experience and lacks quantitative monitoring standards throughout the entire process.

[0010] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application.

[0011] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of this application described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, apparatus, product, or server that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or modules not explicitly listed or inherent to such processes, methods, products, or devices.

[0012] It should be noted that in this specification, "feature point" and "key point" refer to the same object, both representing sparse matchable points used for cross-frame tracking, cross-modal constraint construction, and registration solution.

[0013] Example 1, as Figure 1 As shown, this application provides a real-time analysis method for interventional ultrasound images based on feature point extraction and elastic registration. The method includes: S100: Acquire reference volume data and generate reference safety boundary volume data; It should be noted that the reference volume data can be CT / MRI / 3D ultrasound volume data acquired preoperatively, or 3D ultrasound volume data obtained through free-hand scanning during the intraoperative or preoperative phases. The reference volume data is denoted as... Its coordinate system is The baseline data may further include organ segmentation data. Lesion segmentation data With security boundary data .

[0014] In a preferred implementation, the security boundary volume data Data from lesion segmentation Preset distance for outward expansion And subject to organ segmentation data Restricted generation, among which Preferred selection Specifically, for Data of the expanded lesion were obtained by performing three-dimensional morphological dilation. , and take: This ensures that the safety boundary lies within the organ's perimeter. , and This data serves as input for subsequent rigid / elastic registration and intraoperative monitoring analysis.

[0015] S200: Continuously receive two-dimensional ultrasound image frames output by the interventional ultrasound probe and probe pose data synchronized with the two-dimensional ultrasound image frames in time, perform preprocessing and spatiotemporal calibration on the two-dimensional ultrasound image frames, and generate a spatiotemporally consistent interventional ultrasound image sequence. It should be noted that the interventional ultrasound system in this embodiment includes an ultrasound host, an interventional ultrasound probe, a positioning sensor, and a computing unit. The positioning sensor is an electromagnetic sensor or an optical sensor, which is fixedly installed on the interventional ultrasound probe to output probe pose data. The computing unit is a workstation or embedded computing platform with GPU acceleration capabilities, used to perform image processing, feature point extraction, registration solving, and result display.

[0016] To ensure the feasibility and spatiotemporal consistency of subsequent registration and analysis, the interventional ultrasound image sequence output in this step should meet the following requirements: each frame of two-dimensional ultrasound image has a consistent timestamp with its corresponding probe pose data, and each frame of two-dimensional ultrasound image can obtain a unique spatial pose in the patient coordinate system.

[0017] Specifically, preprocessing and spatiotemporal calibration include the following sub-steps: S201: Perform timestamp alignment and interpolation synchronization on two-dimensional ultrasound image frames and probe pose data to generate one-to-one time-aligned frame pairs; Because the video output link of the ultrasound host and the data output link of the positioning sensor have different sampling frequencies and link delays, there may be a problem that the image frames and pose data do not arrive synchronously. Therefore, the received two-dimensional ultrasound image frames and probe pose data are all appended with original timestamps and written into the image queue and pose queue respectively, and alignment is achieved by timestamp matching.

[0018] In this embodiment, the two-dimensional ultrasound image frame sequence is assumed to be: The probe pose data sequence is as follows: in: For the first Two-dimensional ultrasound images; For the first Timestamp of a frame of two-dimensional ultrasound image; For the first Each probe pose data point must include at least translation and rotation components; For the first Timestamps of probe pose data.

[0019] For each frame of two-dimensional ultrasound image Search the pose queue for the one that satisfies the pose. Adjacent pose data and The pose was obtained by interpolating linear interpolation and spherical linear interpolation. This generates time-aligned frame pairs. To ensure a unique implementation path, the pose queue is stored and maintained in ascending order of timestamps, and is set as a fixed-capacity queue with a maximum length of [value missing]. Preferred selection Each received pose data When the time is up, insert it at the end of the queue according to the timestamp; when the queue length exceeds When the time comes, delete the pose data corresponding to the earliest timestamp at the head of the queue. Assume there are a total of [number] poses in the updated pose queue. The pose data is denoted as ,in and order , . Therefore, when The current time range covered by the pose queue There must exist two adjacent pose data that satisfy certain conditions. To ensure the feasibility of interpolation synchronization, the two adjacent pose data must satisfy... When detected season .when season ;when season ,in for Time-of-flight pose data, for Time-position pose data; and this It forms a time-aligned frame pair with the corresponding two-dimensional ultrasound image frame for output, eliminating the need to discard or wait for pose processing of the two-dimensional ultrasound image frame.

[0020] S201-1: Decompose the probe pose into rotational and translational components. Use linear interpolation for the translational components and spherical linear interpolation for the rotational components. Among them, position It can be represented as a combination of rotation and translation components. ,in For rotational components, use rotation matrices or unit quaternions. Let be a translation component and a three-dimensional translation vector. , And define the interpolation weights: Linear interpolation is then used for the translation components: Spherical linear interpolation (Slerp) is used for the rotational components. When the rotational components are represented as quaternions, let... They are respectively The corresponding unit quaternion (if the input is a rotation matrix, it is first converted to a unit quaternion).

[0021] S201-2: When the dot product of the rotation quaternions of two adjacent pose data is negative, take the opposite of the rotation quaternion of the next pose data to select the shortest rotation arc. Available When Define intermediate quantities And define the included angle: in To preset a minimum threshold, the preferred value is... radians. When Greater than the preset minimum threshold At that time, Slerp is calculated as: S201-3: When the rotation angle is less than the preset threshold, normalized linear interpolation is used for the rotation component.

[0022] when To avoid numerical instability, normalized linear interpolation is used: And After normalization, it is converted into a rotation matrix. The final interpolated pose is obtained as follows: S202: Perform denoising, gain normalization and scale resampling on the two-dimensional ultrasound image frame to generate a standardized two-dimensional ultrasound image frame; Due to speckle noise and gain differences in ultrasound images, denoising and normalization are performed on the two-dimensional ultrasound image frames to improve the stability of feature point extraction. Denoising can employ anisotropic diffusion, nonlocal means, or a convolutional neural network-based denoising network; gain normalization can use quantile stretching or histogram matching to map pixel intensity to a fixed range. .

[0023] To establish a reusable mapping between pixel coordinates and physical coordinates, scale resampling is performed on the two-dimensional ultrasound image frames to unify the pixel spacing to a preset physical scale. ,For example Thus, standardized two-dimensional ultrasound image frames are obtained. .

[0024] Specifically, the raw pixel spacing of a two-dimensional ultrasound image frame is read. ,in , These are the physical distances, one horizontal and one vertical, preferably obtained from imaging parameters or image head information output by the ultrasound equipment. The target pixel spacing is set to... And calculate the scaling factor: Let the original image size be... The size of the resampled image is: in This indicates rounding to the nearest integer. Bilinear interpolation is used for grayscale reconstruction after resampling; when... or First, a two-dimensional Gaussian low-pass filter is applied to the original image to suppress aliasing, where the filter standard deviation is preferably selected as... , (Unit: pixels). Boundary copying and padding are used for sampling locations that exceed the original boundaries. After resampling, a normalized two-dimensional ultrasound image frame is obtained. This unifies the subsequent pixel-to-physical mapping scale, resulting in a pixel-to-physical mapping matrix. The scale term in the and the said Consistent.

[0025] S203: Based on the probe-sensor extrinsic parameter calibration matrix and the sensor-patient coordinate system transformation matrix, calculate the pose prior and generate a spatiotemporally consistent interventional ultrasound image sequence; To map two-dimensional ultrasound image frames to the patient coordinate system, the following coordinate system is defined: Image coordinate system The origin is the top left corner of a standardized two-dimensional ultrasound image frame, and the unit is pixels. Probe coordinate system : The mechanical structure fixed to the interventional ultrasound probe; Sensor coordinate system Fixed to electromagnetic / optical sensors; Patient coordinate system : Fixed to the patient's body surface or operating table.

[0026] In this embodiment, homogeneous coordinates are used to represent three-dimensional points, and the column vector right multiplication convention is adopted. Homogeneous transformation matrix. Indicates the coordinate system Transform the points below to the coordinate system Below, that is, for any point ,have ,in for matrix.

[0027] The extrinsic parameter matrix from the probe coordinate system to the sensor coordinate system is obtained through calibration. The pose matrix from the sensor coordinate system to the patient coordinate system is obtained in real time through the tracking system. Then the first The probe pose matrix of a frame-normalized two-dimensional ultrasound image in the patient coordinate system is: in Timestamps and Image Timestamps The correspondence is determined by the timestamp alignment and interpolation synchronization in step S201, that is, interpolating the sensor pose data at adjacent timestamps to obtain the time. The corresponding sensor pose, and the rotation and translation of the interpolated pose are assembled into a homogeneous matrix representation. .

[0028] Among them, pose prior refers to the current frame timestamp for this step. The obtained probe pose Used for the neighborhood and region of interest of the tip of subsequent instruments. The position and extent provide spatial constraints; the spatial gating involved in subsequent steps can also reuse this pose prior.

[0029] Furthermore, any pixel in a standardized two-dimensional ultrasound image frame Physical coordinates in the image coordinate system A pixel-to-physical mapping matrix get: in By pixel pitch The image origin offset and imaging plane orientation are uniquely determined. Specifically, let the selected physical mapping origin in the pixel coordinate system be... And let its corresponding three-dimensional point in the probe coordinate system be... Let the two unit direction vectors of the imaging plane in the probe coordinate system be respectively and And satisfy , Then any pixel point The three-dimensional physical coordinates in the probe coordinate system are: Accordingly, the mapping matrix can be written as: in , , This is a constant obtained through calibration. In a preferred implementation, the following is selected: The coordinate axes make the imaging plane located at and , Relative to the probe coordinate system axis, If the axes are aligned, then... And the above formula can be degenerated into , .

[0030] Finally, the three-dimensional position of this pixel in the patient coordinate system is: Through the above-mentioned time synchronization, image standardization, and spatial pose calculation, a spatiotemporally consistent interventional ultrasound image sequence is generated: Each frame contains a standardized two-dimensional ultrasound image and the corresponding patient's pose in the coordinate system.

[0031] S300: Continuously read the interventional ultrasound image sequence, extract the feature points of anatomical structures and interventional instruments, construct feature point map information, and generate a feature point map information sequence with confidence. It should be noted that ultrasound images may be obscured by acoustic shadowing, artifacts, and high echoes from instruments during intervention. Relying solely on pixel-level intensity similarity can easily lead to registration drift. Therefore, this embodiment extracts feature points from both anatomical structures and interventional instruments simultaneously and encodes the topological relationships between these feature points into feature point map information to enhance the stability and interpretability of the matching.

[0032] In this embodiment, the feature point map information sequence output should include at least: key point coordinates, key point confidence, key point descriptor, key point category (anatomical structure / instrument), and edge set and geometric consistency constraints constructed based on the key points.

[0033] Specifically, feature point extraction and feature point map construction include the following sub-steps: S301: Determine the region of interest based on the segmentation mask and pose prior, and generate a set of image patches for the region of interest; Standardized two-dimensional ultrasound image frames Perform organ contour segmentation and instrument segmentation to obtain organ segmentation masks. Instrument segmentation mask In one implementation, the segmentation model is applied to the input... Output organ probability map Device Probability Chart and using a fixed threshold Binarization: in Preferred selection To suppress artifacts and isolated noise, a morphological closing operation is performed on the binary mask, and areas smaller than a threshold are removed. Connected components, Preferred selection Pixels; the connected component with the largest preserved area in the organ mask is used as the organ's outer contour region. To ensure cross-frame stability, the execution time consistency between the current frame mask and the previous frame mask is constrained, and the radius of the structuring element is set to... The dilation operator is ,but: in Preferred selection Pixel; when season and If the above update results in the mask being an empty set, then respectively set the mask to an empty set. , ,when And when it is still an empty set, let them be respectively and It is an empty set.

[0034] The region of interest is determined based on the segmentation mask. ,in It should at least cover the organ boundary and the vicinity of the instrument tip. Specifically, the boundary zone region of the organ mask should be calculated. ,in Preferred selection Pixel; calculate the skeleton lines of the machine mask and take the skeleton endpoints as candidate tip points, then select the endpoint with the smallest Euclidean distance to the tip point of the previous frame as the tip point of the current frame. When there are no skeleton endpoints, let The centroid of the instrument mask. Construct a structure with a side length of centered on . square neighborhood , Preferred selection Pixels. Let the joint region be: and take The axis-aligned minimum bounding rectangle is used as , and will Cropping to within the effective boundaries of the image; when When the empty set is in order This refers to the entire image region. Cropped into a set of image patches of a fixed size: The image patch size is fixed as The step size is fixed. , Preferred selection , Preferred selection ;when When the boundary is insufficient to crop a complete image patch, boundary copying and padding are used, and This represents the number of image blocks.

[0035] S302: Perform keypoint detection on the set of image patches in the region of interest to generate a set of feature points; For each image patch The keypoint detection process outputs a keypoint heatmap, a keypoint offset map, and a descriptive subgraph, which are then further decoded to obtain a set of keypoints. in: For the key point Two-dimensional coordinates in; For key point confidence, ; Keypoint descriptor subvectors are used for cross-frame matching; These are key point category labels used to distinguish between anatomical key points and instrument key points.

[0036] S302-1: Output keypoint heatmap, keypoint offset map and descriptor map for each image block; To ensure the uniqueness of the keypoint detection and decoding path, the keypoint detection model uses image patches. As input, output three types of feature maps: keypoint heatmaps. Key point offset map With descriptor graph .in , , And the network downsampling step size is fixed at 1. ,Right now , The descriptor dimension is fixed as .

[0037] S302-2: Perform local maximum suppression on the key point heatmap to determine the position of candidate key points, and perform sub-pixel decoding on the position of candidate key points based on the key point offset map to obtain the key point coordinates; heat map Perform local maximum suppression and retain the values ​​that satisfy the condition. And for Candidate key points for local maxima , where confidence threshold Preferred selection For each candidate keypoint location, in As the confidence level of the key points, the sub-pixel coordinates of the key points in the image patch coordinate system are decoded as follows: The key point coordinates in the image patch coordinate system are transferred through the image patch onto the standardized image. offset of the top left corner Mapped to full-map coordinates: S302-3: Bilinear sampling and normalization are performed on the candidate key point positions on the descriptor graph to obtain key point descriptors; For keypoint descriptors, bilinear sampling is used. Get it from and conduct Normalization: S302-4: The key points obtained across image blocks are deduplicated according to the deduplication distance threshold and filtered according to the confidence level to obtain the set of feature points.

[0038] Key point category labels The unique value is assigned based on the segmentation mask obtained in step S301: when Falling into the instrument segmentation mask Inner Season Otherwise After merging the keypoints of all image patches, select the pairs of points in the overall image whose distance is less than a threshold. The key point is to perform deduplication and retain those with higher confidence, among which Preferred selection Pixels; finally, sort by confidence level from highest to lowest and keep only the top [pixels]. One key point, Preferred selection Thus, the set of feature points is obtained. .

[0039] S303: Perform a matching score calculation on the feature point set to generate a candidate matching pair set; To perform cross-frame tracking and cross-modal constraint construction, for adjacent frames Or, for the current frame feature point set and the reference volume data projection feature point set. Calculate the descriptor distance and form a set of candidate matching pairs: in The matching score is used to indicate the reliability of the match; a higher score means a more credible match.

[0040] To ensure the uniqueness of the matching score calculation path, the two sets of feature points to be matched are denoted as follows: and And all of the descriptors have been... Normalization. Candidate matches are computed only when the keypoint categories are consistent, i.e., when... Descriptor similarity is calculated only when the categories are inconsistent; when the categories are inconsistent, descriptor similarity is calculated only when the categories are inconsistent. For keypoint pairs with consistent categories, define cosine similarity: And only retain those that meet the requirements. Candidate pairs, where the threshold Preferred selection Furthermore, for each Among the candidates that meet the threshold, the first nearest neighbor with the highest similarity is selected. With the second nearest neighbor And calculate the ratio: in To prevent extremely small constants with a denominator of zero, it is preferable to take... ;when Keep matching pairs ,in Preferred selection To improve the uniqueness and interpretability of the matching, a nearest neighbor constraint is further adopted: only if... and Only when they are established at the same time will Added as a candidate matching pair .

[0041] To suppress unreasonable matching by utilizing pose priors, a pixel displacement threshold is set for matching scenes between adjacent frames. For a fixed value, the preferred value is... Pixels; only when The candidate match is retained. For cross-modal matching scenarios, the reference volume data is projected onto the feature point set. Projecting onto the current frame plane yields reference two-dimensional coordinates. And the same threshold is used for spatial gating. Finally, for each candidate matching pair The overall matching score is defined as follows: And all candidate matches by Sort from highest to lowest, keeping only the top few. There are 10 candidate matching pairs, of which Preferred selection .

[0042] S304: Perform outlier removal based on feature point map geometric consistency constraints to generate a set of corrected matching pairs; To reduce registration drift caused by mismatches, an undirected graph is constructed based on the feature point set: in: A set of nodes; The graph is a set of edges, connecting pairs of nearby key points in space, and geometric properties are computed for each edge. To ensure the uniqueness of the graph construction path, for each node... Within the same set of key points (i.e.) Selecting its Euclidean distance as the metric The nearest neighbors form an undirected edge. Preferred selection And only retain those with side lengths not exceeding the threshold. The edge, Preferred selection Pixels. For any edge Define the edge vector and the edge length as The direction angle of the side is Simultaneously, to calculate the local texture gradient direction difference, the image gradient direction at the keypoint is taken as... And define the gradient direction difference as ,in Standardize the angle to .

[0043] For the set of candidate matching pairs Outlier removal is performed using the geometric consistency constraint of the edge set: when the rate of change of the matched edge length exceeds a threshold. Or the difference in side direction exceeds the threshold The matching pair is identified as an outlier and removed, thus obtaining a set of corrected matching pairs. .

[0044] Specifically, let the set of candidate matching pairs be... ,in For the index of key points in the current frame, For reference to the key point index, and for each Keep at most one For any matching pair Statistics on its position in the graph Neighborhood consistency in: Let For nodes The set of neighbors, and only considers neighbor nodes that both exist simultaneously. Its corresponding match is For each such neighbor pair, calculate the rate of change of length and direction difference between the current frame edge and the reference edge: in Indicates reference key points Coordinates in the same two-dimensional coordinate system (coordinates of the previous frame when matching adjacent frames, and coordinates projected onto the current frame plane when matching across modalities). It is a minimal constant, and preferably takes .when and When a pair of neighbors is matched, that neighbor pair is counted as a match. If the number of its consistent neighbors is less than the threshold These are then identified as outliers and removed. Preferred selection , Preferred selection , Preferred selection After removal, the remaining matching pairs are scored according to their matching characteristics. Sort from highest to lowest, keeping only the top few. Each matching pair is considered as a set of corrected matching pairs. ,in Preferred selection .

[0045] The final output is a sequence of feature point map information with confidence scores: S400: Continuously read the feature point map information sequence and the reference volume data, perform rigid registration and elastic registration on the reference volume data based on the feature point map information, generate a rigid transformation matrix and a three-dimensional elastic deformation field, and generate the deformed reference volume data and the deformed safety boundary. It should be noted that the reference body data The corresponding lesion segment data and safety boundary data have been acquired and generated in step S100, and are used as input data for registration solution in this step.

[0046] To meet real-time requirements, this embodiment adopts a combined strategy of "initial rigid alignment based on pose prior + elastic compensation driven by feature point map constraints + periodic secondary fine registration". Elastic deformation is modeled through a three-dimensional elastic deformation field, so that the reference body data dynamically fits the deformation state during surgery in the patient coordinate system.

[0047] Specifically, rigid registration and flexible registration include the following sub-steps: S401: Perform an initial rigid transformation on the reference body data based on the probe pose data to generate initial aligned body data; Let the reference body data be Its coordinate system is Obtained through preoperative registration or intraoperative initialization. To the patient coordinate system The initial rigid transformation matrix , and combined with the first Frame probe pose matrix It can generate spatial relationships for extracting reference slices.

[0048] in for The homogeneous rigid transformation matrix is ​​used to transform points in the reference coordinate system. Transform to the lower point of the patient coordinate system ,satisfy: In a preferred implementation, the Through rigid registration of point sets, at least [number] points were collected in the patient coordinate system. anatomical landmarks or surface landmarks And extract the corresponding marker points in the reference body data coordinate system. ,in And preferred By minimizing Solving for rotation and translation yields .

[0049] For a consistent description, the rigid transformation matrix is ​​represented as: in For rotational components, This represents the translation component.

[0050] Perform an initial rigidity transformation on the reference volume data to obtain the initial aligned volume data: in Take as Used to transfer reference body data from Transform to In order to make The implementation path is uniquely determined, and the sampling grid of the output volume data is defined in the patient coordinate system. and to Center point of any voxel Calculate its corresponding point in the reference coordinate system: Again exist Trilinear interpolation sampling was performed at the location to obtain ;when Falling When the value is outside the effective range, a constant value is used for filling, and the filling value is preferably selected. .

[0051] To generate a reference slice corresponding to the real-time two-dimensional ultrasound image, the transformation matrix from the probe coordinate system to the reference coordinate system is calculated: Based on this, the two-dimensional imaging plane determined by the pixel-to-physical mapping relationship described in step S203 is mapped onto the reference volume data or the initial alignment volume data to extract the reference slice.

[0052] S402: Construct a registration objective function based on feature point map information to generate an optimizable joint constraint expression; To simultaneously utilize the consistency of sparse feature points and the consistency of image gray levels, a registration objective function is constructed: in: This is the feature point reprojection error term, used to constrain the error in projecting the deformed reference feature points in the reference volume data onto the current two-dimensional ultrasound image; This is an image similarity term used to constrain the similarity between simulated slices of deformed volume data on the current imaging plane and real two-dimensional ultrasound images; This is a deformation smoothing term used to constrain the gradient change of the elastic deformation field; This is a locally approximate rigid constraint used to keep the relative distance changes in the local neighborhood within a reasonable range. The weighting coefficients satisfy... and and select the best ; In the reference volume data, the reference feature point refers to the feature point detected on the reference slice extracted based on the current frame probe pose in the reference volume data (or the initial aligned volume data obtained in step S401). After mapping the reference feature point on the reference slice to three-dimensional space according to the pixel-to-physical mapping and spatial pose relationship described in step S203, its three-dimensional coordinates in the reference coordinate system are obtained, denoted as . The correspondence between reference feature points and key points in the current two-dimensional ultrasound image is determined by the corrected set of matching pairs output from steps S303 and S304. Provided.

[0053] In this embodiment, the feature point reprojection error term can be defined as: in: These are the coordinates of key points in the current two-dimensional ultrasound image; The three-dimensional coordinates of the reference key points in the reference volume data; This represents the coordinate mapping after the application of a three-dimensional elastic deformation field. This represents the projection operator from three-dimensional coordinates to the current two-dimensional imaging plane; The confidence level of key points is used to weight and suppress the influence of low-confidence key points.

[0054] To enable the projection operator The implementation path is uniquely determined, making Let be the rigid transformation matrix from the reference coordinate system to the patient coordinate system, and Indicates reference key points Coordinates in the patient coordinate system. For any three-dimensional point in the patient coordinate system. First, calculate its coordinates in the probe coordinate system: Then use the pixel-to-physical mapping parameters defined in step S203. , , , With pixel spacing Calculate its pixel coordinates in the image coordinate system: in for The first three components.

[0055] The image similarity term is defined using normalized cross-correlation (NCC): in This means extracting the corresponding slice from the volume data according to the current probe pose and generating a simulated two-dimensional image. To make... The implementation path is uniquely determined, so that the output simulated slice and With the same resolution and pixel pitch, for any pixel point According to step S203 and Its three-dimensional position in the patient's coordinate system can be obtained: Then analyze the volume data after deformation. exist The original simulated grayscale is obtained by performing trilinear interpolation sampling at the location. ;when When it falls outside the valid range of volume data, let .right Performing linear normalization yields : in , , It is a minimal constant, and preferably takes and order Furthermore, The region of interest obtained in step S301 pixel set within Let the computational domain be defined as: in and They are respectively and exist The mean of the above.

[0056] S403: Perform elastic deformation field estimation based on deformation parameterization to generate a three-dimensional elastic deformation field; S403-1: Parameterizing the deformation field and determining the effective deformation domain based on cubic B-spline control point grid; To balance real-time performance and controllability, the three-dimensional elastic deformation field is parameterized as a B-spline control point grid. The control point grid is defined in the patient coordinate system. Below, a cubic B-spline basis function is used for construction. Let the control point grid be... Each control point has a three-dimensional displacement vector. Then any three-dimensional point The transformation is: in These are B-spline basis functions.

[0057] To ensure the control point grid construction path is uniquely determined, let the effective deformation domain be... It is taken as the region of interest in the current frame. The bounding box and its outer radius in the patient coordinate system The union of, where Preferred selection .exist Construct an equally spaced grid of control points, with a spacing of [missing information]. Preferred selection Furthermore, the control point grid is expanded outward by one control point spacing in each of the three axes to ensure that deformation at the boundary is definable. The control point positions can be represented as: in Integer index, This is the starting control point after expanding the bounding box from its smallest corner point. For any point... Only for its own grid cell Sum of adjacent control points, and It is uniquely determined by the tensor product of the one-dimensional basis functions of the cubic B-spline.

[0058] S403-2: Apply discrete bending energy smoothing constraints to control point displacements to suppress discontinuous deformation; The smoothing term is defined as a discrete bending energy constraint on the displacement of the control points: The second-order difference operator is defined as follows: , and Similarly, the summation range is the index of all control points with valid second-order difference neighborhoods.

[0059] S403-3: Constructing local approximate rigid constraints based on the K nearest neighbors of reference key points to maintain neighborhood distance consistency; The local approximate rigid constraint is defined on the set of reference key points used for registration. Let the set of reference key points be... And it is obtained by rigid transformation and deformation mapping onto the patient coordinate system. For each Using Euclidean distance as a metric in a set Select it from within The nearest neighbors are used as the neighborhood. ,in Preferred selection and only retain distances not exceeding Neighboring points, Preferred selection A local approximate rigid constraint can be defined as: Among them when If the set is empty, skip this item.

[0060] Through the Perform iterative optimization to obtain the optimal result. With control point displacement set This generates a three-dimensional elastic deformation field. To ensure the optimization process is feasible and the path is uniquely determined, the initialization command is set to... and order When there is a deformation field from the previous frame. At that time, the values ​​are sampled at equal intervals to the current control point grid as initial values. Unconstrained optimization based on L-BFGS is used for updating. The maximum number of iterations is Preferred selection When satisfied Stop iteration when, where Preferred selection To suppress unreasonable large deformations, amplitude clipping is applied to the displacement of each control point. ,in Preferred selection After trimming, the iteration continues. After the iteration ends, the displacement of the control points is calculated using the same B-spline interpolation to obtain the displacement vector at any voxel, thus obtaining the three-dimensional elastic deformation field defined in the patient coordinate system. .

[0061] S404: Apply a three-dimensional elastic deformation field to the reference body data and safety boundary to generate deformed reference body data and deformed safety boundary. Applying a three-dimensional elastic deformation field to the initial aligned body data yields the deformed reference body data: The same deformation field was applied to both the lesion volume data and the safety boundary volume data to obtain the deformed lesion volume data. With the deformed safety boundary volume data .

[0062] In order to make The implementation path is uniquely determined, and a reverse mapping resampling strategy is adopted. Let Sampling grid defined in patient coordinate system Above, and Represents the lower point of the patient's coordinate system The three-dimensional displacement vector (unit: For any output voxel center point Calculate the sampling points corresponding to it before deformation: And on exist Trilinear interpolation sampling was performed at the location to obtain .when Exceeding When the effective range is reached, mirror boundary filling is used. Interpolation sampling is performed after each coordinate component is mapped back to the valid range.

[0063] For both lesion region volume data and safety boundary volume data, a binary mask representation is preferred. To maintain the boundary stability of discrete labels, the same inverse mapping is used... and Resampling, and using nearest neighbor interpolation for the mask: for any ,exist Take the nearest neighbor voxel value and set it to be greater than or equal to 1. Set the value to The rest are set to The deformed binary mask is obtained. and ;when When the boundary exceeds the valid range, mirror boundary filling is also used.

[0064] S405: Set the cache and trigger secondary fine registration to generate a stable and updated registration result; To reduce the impact of single-frame noise on registration, the queue attribute is configured to store the most recently used frames. Frame data, where the buffer length For a fixed positive integer, the preferred choice is... : To ensure that the triggering and update paths are uniquely determined, the triggering period for the secondary fine-registration is set to a fixed time interval. Preferred selection When the trigger time is reached and the number of valid frames in the buffer reaches 100,000, the trigger time is reached. At that time, a second fine-registration update is performed. The second fine-registration is based on the estimate of the current frame. As initialization, the registration loss of each frame in the buffer sequence is accumulated as the optimization objective: in , , and The calculations are performed according to the definitions of steps S402 and S403, respectively, but the input of the current frame is replaced with the input of the first frame in the buffer. Frames , and The same optimizer and stopping conditions as in step S403 are used. Perform iterative updates, and after the iteration is complete, the updated version will be... This serves as the initial value for the next frame, thereby improving stability and suppressing cumulative drift.

[0065] S500: Continuously reads real-time two-dimensional ultrasound images and deformed reference body data, performs interventional device positioning and lesion ablation status analysis, generates monitoring indicators and risk warning information and outputs them; It should be noted that this step uses the deformed reference body data output in step S404. With the deformed safety boundary volume data As a "reference after intraoperative deformation", it generates intuitive and quantifiable monitoring indicators on the two-dimensional ultrasound real-time image, including the three-dimensional position of the instrument tip, the distance from the tip to the safety boundary, the ablation coverage rate, and the ablation completion prompt.

[0066] Specifically, the analysis of interventional device localization and lesion ablation status includes the following sub-steps: S501: Based on the instrument segmentation mask and instrument feature points, perform instrument axis fitting to generate instrument tip position and three-dimensional positioning results; Within the two-dimensional ultrasound image plane, based on instrument segmentation mask The skeleton line of the hyperechoic connected region of the instrument is extracted, and combined with key points of the instrument category in the instrument feature point set, robust straight line fitting is performed on the skeleton line to obtain the instrument axis. With two-dimensional tip To ensure a unique implementation path, let the standardized two-dimensional ultrasound image size be [value missing]. and adopt Neighborhood defines connectivity. (For) Perform connected component analysis and retain the connected component with the largest area. ;when When it is an empty set, let Take the tip of the previous frame ,when season Center point of the image And skip subsequent axis fitting.

[0067] right Execution refinement yields the skeleton mask. And extract the skeleton pixel set. Define the skeleton endpoint set. .when When not empty, calculate the distance from each endpoint to the image boundary. and the endpoint The smallest is defined as the instrument base endpoint. , The largest is defined as the initial value at the instrument tip. ;when When it is an empty set, let and Take respectively The centroid and the skeleton point furthest from the centroid.

[0068] To fit the instrument axis, a two-dimensional point set is constructed for fitting. It is a set of skeleton points. The union of the coordinates of key points of category "device" in the set of device feature points. Let the set of device feature points be... And assign weights to each instrument feature point. Assign a fixed weight to each skeleton point. Preferred selection Calculate the weighted mean: And calculate the weighted covariance matrix: Pick The unit eigenvector corresponding to the largest eigenvalue As the direction of the axis. and Determine the direction sign, if Then let The axis of the instrument is represented as: The two-dimensional apex is defined as the orthogonal projection of the initial apex value onto the axis: Furthermore, in conjunction with the first The frame probe pose matrix and pixel-to-physical mapping matrix map the two-dimensional tip point to the three-dimensional tip position in the patient coordinate system: S502: Performs ablation region segmentation based on real-time two-dimensional ultrasound images and generates an ablation region mask; Standardized two-dimensional ultrasound image frames Perform ablation region segmentation to obtain the ablation region mask. To ensure a unique path is determined, a binary classification segmentation model is used to segment the input. Output ablation probability map And only in the region of interest obtained in step S301 Internal calculation and update of the ablation mask; Pixels other than Using a fixed threshold Binarization: in Preferred selection To suppress the porosity and isolated points caused by speckle noise, [the following measures were taken]. Perform morphological closing and opening operations sequentially, and delete items with areas smaller than the threshold. Connected components, Preferred selection Pixel.

[0069] S503: Project the deformed safety boundary onto the two-dimensional ultrasound image plane and calculate the ablation coverage rate to generate an ablation completion determination result; The deformed safety boundary body data The projection safety boundary mask is obtained by projecting the current probe pose onto a two-dimensional plane. To ensure the projection calculation path is uniquely determined, only the region of interest obtained in step S301 is considered. Endogenous generation and will Outside pixels set to For any pixel According to step S203 and Its three-dimensional position in the patient's coordinate system can be obtained: Then, the binary safe boundary volume data exist Nearest neighbor interpolation sampling and binarization are performed at the location to obtain: when Exceeding When the effective range is reached, let Ablation coverage is defined as: in Indicates to The area in pixels is obtained by summing the pixel values ​​of the inner binary mask. To prevent extremely small constants with a denominator of zero, it is preferable to take... ;when season .

[0070] To ensure that the completion decision path is uniquely determined, a coverage threshold is set. Stable frame rate threshold Define a stable counter. ,when season Otherwise ,in ;when When the ablation is complete, a notification message will be output.

[0071] S504: Calculates risk warning information based on the distance from the tip to the safety boundary and integrates and displays the output; Based on three-dimensional tip position With the deformed safety boundary volume data Calculate the minimum distance from the tip to the safety boundary: in For the first The instrument tip at frame time in the patient coordinate system The three-dimensional position vector below; For the first The deformed safety boundary volume data corresponding to the frame is preferably a binary voxel mask defined in the patient coordinate system; This represents the set of boundaries of the security boundary volume data; For any boundary point in the boundary set; This is the minimum Euclidean distance from the instrument tip to the safety boundary; It is a three-dimensional Euclidean norm.

[0072] To ensure that distance calculation and path determination are uniquely determined, Represented as a binary mask for the safety boundary after deformation. The boundary voxel set, that is, satisfying And in There is at least one voxel in the neighborhood. Make The set of voxel centers. For Construct a point set from the center points of all boundary voxels. And it is calculated using three-dimensional Euclidean distance: in and The coordinates of the voxel center point in the patient coordinate system; Neighborhood representation on a 3D voxel mesh A set of voxels adjacent in a single axial direction; It is a three-dimensional point set composed of the center points of boundary voxels; the For point set The three-dimensional coordinates of any point in the array, and with respect to... They are placed in the same coordinate system to ensure that the dimensions of distance calculations are consistent.

[0073] To reduce computational load, it is preferable to construct... When performing voxel spacing is Equal-interval downsampling, Preferred selection To ensure the risk assessment path is uniquely determined, a distance threshold is set. ;when A risk warning message is output when the warning is triggered; otherwise, no risk warning message is output. The sampling interval (in units) for downsampling the boundary point set. ); Distance threshold for risk warning (unit: ).

[0074] Output methods include: fusing and displaying corresponding slices of real-time two-dimensional ultrasound images with deformed reference body data, and overlaying lesion boundaries, safety boundaries, instrument axis, instrument tip, ablation coverage values, and risk warning information.

[0075] S600: Continuously monitors registration confidence and abnormal states, and performs degradation and re-initialization control when abnormal states are identified; It should be noted that during the intervention process, situations may occur such as the probe leaving the target area, strong instrument obstruction, or a sudden drop in image quality, which may reduce the number of feature points or make the matching unreliable. Therefore, it is necessary to conduct a confidence assessment of the registration results and provide feasible degradation strategies.

[0076] In this embodiment, the registration confidence index is defined as follows: in The mean confidence score of keypoints in the corrected matching pair set. As weight.

[0077] To ensure that the triggering conditions for downgrades and reinitializations are uniquely determined, a confidence threshold is set. Minimum number of valid matching pairs threshold And reinitialize the threshold for consecutive low-confidence frames. ,in The selectable range is 0.4-0.8, preferably 0.5. The selectable range is 50–150, preferably 80. The selectable range is 20–60 frames, preferably 30 frames. When Or by adjusting the number of matched pairs When this occurs, a degradation strategy is implemented: the elastic deformation field update is frozen, and only the rigid transformation matrix update is retained; when the conditions are continuously met... or The frame rate reached When this happens, the re-initialization process is triggered, and the initial alignment of the reference volume data to the patient coordinate system is re-executed and the feature point map information is re-established.

[0078] During use, the operator activates the interventional ultrasound real-time analysis system. The system continuously receives two-dimensional ultrasound image frames and probe pose data, performs preprocessing and spatiotemporal calibration on the two-dimensional ultrasound image frames to generate a spatiotemporally consistent interventional ultrasound image sequence. The system extracts anatomical structures and interventional instrument feature points from each image frame and constructs feature point map information. Based on the feature point map information, it performs rigid / elastic joint registration on the reference body data to generate deformed reference body data and safety boundaries. The system further outputs instrument positioning results, distance from the tip to the safety boundary, ablation coverage, and risk warning information on the real-time two-dimensional ultrasound screen, and executes degradation and reinitialization strategies when registration is unstable, thereby achieving real-time stable monitoring and quantitative evaluation of the entire interventional procedure.

[0079] It should be noted that the order of the embodiments described above is merely for descriptive purposes and does not represent the superiority or inferiority of the embodiments. Furthermore, the above description focuses on specific embodiments of this specification. Additionally, the processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired results. In some implementations, multitasking and parallel processing are possible or may be advantageous.

[0080] The above description is only a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

[0081] This specification and accompanying drawings are merely illustrative examples of this application and are intended to cover any and all modifications, variations, combinations, or equivalents within the scope of this application. Clearly, those skilled in the art can make various alterations and modifications to this application without departing from its scope. Therefore, if such modifications and variations fall within the scope of this application and its equivalents, this application intends to include such modifications and variations.

Claims

1. A real-time analysis method for interventional ultrasound images based on feature point extraction and elastic registration, characterized in that, The method includes: Acquire reference volume data and generate reference safety boundary volume data; The system continuously receives two-dimensional ultrasound image frames output by the interventional ultrasound probe and probe pose data synchronized with the two-dimensional ultrasound image frames in time. It performs preprocessing and spatiotemporal calibration on the two-dimensional ultrasound image frames to generate a spatiotemporally consistent interventional ultrasound image sequence. The interventional ultrasound image sequence is continuously read, feature points of anatomical structures and interventional instruments are extracted, feature point map information is constructed, and a feature point map information sequence with confidence is generated. The feature point map information sequence and reference volume data are continuously read. Based on the feature point map information, rigid registration and elastic registration are performed on the reference volume data to generate a rigid transformation matrix and a three-dimensional elastic deformation field, and to generate the deformed reference volume data and the deformed safety boundary. Continuously read real-time two-dimensional ultrasound images and deformed reference body data, perform interventional device positioning and lesion ablation status analysis, generate monitoring indicators and risk warning information and output them; Continuously monitor registration confidence and abnormal states, and perform degradation and re-initialization control when abnormal states are identified.

2. The method as described in claim 1, characterized in that, The process of continuously receiving two-dimensional ultrasound image frames output by the interventional ultrasound probe and probe pose data synchronized with the time of the two-dimensional ultrasound image frames, performing preprocessing and spatiotemporal calibration on the two-dimensional ultrasound image frames, and generating a spatiotemporally consistent interventional ultrasound image sequence includes: Perform timestamp alignment and interpolation synchronization on two-dimensional ultrasound image frames and probe pose data to generate one-to-one time-aligned frame pairs. Denoising, gain normalization, and scale resampling are performed on two-dimensional ultrasound image frames to generate standardized two-dimensional ultrasound image frames. Based on the probe-sensor extrinsic parameter calibration matrix and the sensor-patient coordinate system transformation matrix, the pose prior is calculated and a spatiotemporally consistent interventional ultrasound image sequence is generated.

3. The method as described in claim 2, characterized in that, The continuous reading of the interventional ultrasound image sequence, extraction of feature points of anatomical structures and interventional instruments, construction of feature point map information, and generation of a feature point map information sequence with confidence include: Regions of interest are determined based on segmentation masks and pose priors, and a set of image patches of regions of interest is generated. Perform keypoint detection on a set of image patches in the region of interest to generate a set of feature points; Perform a matching score calculation on the feature point set to generate a candidate matching pair set; Outlier removal is performed based on the geometric consistency constraints of the feature point map, generating a set of corrected matching pairs.

4. The method as described in claim 1, characterized in that, The continuous reading of the feature point map information sequence and the reference volume data, the rigid registration and elastic registration of the reference volume data based on the feature point map information, the generation of a rigid transformation matrix and a three-dimensional elastic deformation field, and the generation of deformed reference volume data and deformed safety boundaries include: Based on the probe pose data, an initial rigidity transformation is performed on the reference body data to generate initial aligned body data; Based on feature point map information, a registration objective function is constructed to generate an optimizable joint constraint expression; Elastic deformation field estimation is performed based on deformation parameterization to generate a three-dimensional elastic deformation field. A three-dimensional elastic deformation field is applied to the reference volume data and the safety boundary to generate the deformed reference volume data and the deformed safety boundary. Set up a cache and trigger a second fine registration to generate a stable and updated registration result.

5. The method as described in claim 1, characterized in that, The system continuously reads real-time two-dimensional ultrasound images and deformed reference body data, performs interventional device positioning and lesion ablation status analysis, generates monitoring indicators and risk warning information, and outputs the following: Based on the instrument segmentation mask and instrument feature points, the instrument axis is fitted to generate the instrument tip position and three-dimensional positioning results; Ablation region segmentation is performed based on real-time two-dimensional ultrasound images to generate an ablation region mask; The deformed safety boundary is projected onto the two-dimensional ultrasound image plane and the ablation coverage is calculated to generate a determination result of ablation completion. Risk warning information is calculated based on the distance from the tip to the safety boundary and then displayed in a fused manner.

6. The method as described in claim 2, characterized in that, The process of performing timestamp alignment and interpolation synchronization between two-dimensional ultrasound image frames and probe pose data includes: The probe pose is decomposed into rotational and translational components. Linear interpolation is used for the translational components, and spherical linear interpolation is used for the rotational components. When the dot product of the rotation quaternion of two adjacent pose data is negative, take the opposite number of the rotation quaternion of the next pose data to select the shortest rotation arc. When the rotation angle is less than a preset threshold, normalized linear interpolation is applied to the rotation component.

7. The method as described in claim 3, characterized in that, The keypoint detection performed on the set of image patches in the region of interest includes: For each image block, output a keypoint heatmap, a keypoint offset map, and a descriptor map; Local maximum suppression is performed on the key point heatmap to determine the location of candidate key points, and sub-pixel decoding is performed on the location of candidate key points based on the key point offset map to obtain the key point coordinates; The candidate keypoint locations are bilinearly sampled and normalized on the descriptor graph to obtain keypoint descriptors; The key points obtained across image blocks are deduplicated according to a deduplication distance threshold and then filtered according to confidence level to obtain the set of feature points.

8. The method as described in claim 4, characterized in that, The process of performing elastic deformation field estimation based on deformation parameterization to generate a three-dimensional elastic deformation field includes: The deformation field is parameterized and the effective deformation domain is determined based on cubic B-spline control point grid. Discrete bending energy smoothing constraints are applied to the displacement of control points to suppress discontinuous deformation; Local approximate rigid constraints are constructed based on the reference key point K nearest neighbors to maintain the consistency of neighborhood distances.

Images