Edge feature point extraction and registration method for rigid body registration of three-dimensional volume data

Abstract

The invention relates to the technical field of medical image processing, and discloses an edge feature point extraction and registration method for rigid registration of three-dimensional volume data, which comprises the following steps: acquiring first and second three-dimensional volume data to be registered; synchronously outputting an edge feature point thermodynamic diagram and a tangent plane attitude parameter diagram by using a full convolutional neural network; analyzing sub-voxel positions of the feature points based on the thermodynamic diagram, and constructing tangent plane normal vectors in a physical space by using attitude parameters; establishing a matching relationship according to the category of the feature points, and constructing a target function based on the distance from the point to the tangent plane; and solving rigid body transformation parameters by adopting a Gaussian-Newton iterative algorithm based on Lie algebra parameterization. According to the method, the geometric constraint mechanism from the point to the tangent plane is introduced, the feature points are allowed to slide in the target tangent plane, the sensitivity of the registration result to the feature point extraction position error is effectively reduced, and therefore the precision and stability of medical image rigid body registration are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, specifically to a method for edge feature point extraction and registration for rigid body registration of three-dimensional volume data. Background Technology

[0002] Medical image registration is a fundamental and crucial technology in medical image analysis and clinical auxiliary diagnosis and treatment. It aims to transform two or more images of the same patient acquired at different times, modally, or from different perspectives into a unified coordinate system. For three-dimensional volumetric data such as CT and MRI, rigid body registration primarily addresses the spatial alignment of anatomical structures with relatively small deformations, such as bones and the brain. Its core task is to solve for the rotation matrix and translation vector that describe the spatial positional differences between the source and target images. Accurate registration results provide physicians with multi-source fusion spatial information, which is of significant clinical importance for developing radiotherapy plans, performing surgical navigation, and monitoring lesion development.

[0003] Existing feature-based 3D rigid body registration techniques typically employ a "feature extraction – matching – transformation solution" process. Specifically, the algorithm first uses edge detection operators or deep learning models to extract representative anatomical landmarks or surface contour point sets from the 3D volume data. After establishing the correspondence between feature points in the source data and those in the target data, an optimization objective function is constructed and solved using the iterative nearest point algorithm or its variants. These methods generally use Euclidean distance as a metric, estimating rigid body transformation parameters by minimizing the coordinate distance between matching point pairs, thereby ensuring that the feature point set of the source image spatially overlaps as much as possible with the feature point set of the target image.

[0004] However, the point-to-point distance constraint mechanism commonly used in existing technologies has limitations in practical applications. Its main problem lies in its excessive sensitivity to the accuracy of feature point extraction locations. When processing medical images affected by partial volume effects or imaging noise, the feature points extracted by the algorithm often struggle to be precisely located at completely identical anatomical positions, inevitably resulting in slight positional deviations between corresponding points in the source and target images. In such cases, strict point-to-point constraints force the source and target feature points to perfectly coincide in coordinates. This mandatory constraint causes the positioning error from the feature extraction stage to be directly introduced into the transformation parameter calculation process, resulting in the final calculated rigid body transformation parameters deviating from their true values, thus reducing the overall registration accuracy. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method for edge feature point extraction and registration for rigid body registration of three-dimensional volume data, which solves the problems mentioned in the background section.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] The first aspect of the present invention provides a method for edge feature point extraction and registration for rigid body registration of three-dimensional volume data.

[0008] This method processes 3D medical volume data stored in voxel matrix form and containing spatial positioning metadata. This invention uses a deep learning network to simultaneously extract the positional and geometric pose information of anatomical edge feature points, and utilizes geometric constraints from points to tangent planes to construct an optimization model, achieving high-precision rigid body registration of source and target data in the physical coordinate system. The method includes the following main processes:

[0009] First, the first 3D volume data to be registered is acquired as the source data, and the second 3D volume data is acquired as the target data. Metadata for both needs to be read to determine their spatial position in a unified coordinate system.

[0010] Secondly, a fully convolutional neural network structure is used to extract features from the input 3D volume data. This network is configured with a dual-branch output structure: the first output branch generates an edge feature point heatmap to characterize the probability distribution of the presence of specific anatomical landmarks in the image space; the second output branch generates a tangential plane attitude parameter map to characterize the geometric orientation of the local surface at the landmark. The attitude parameters consist of angular components in the Euler angle system, including pitch, yaw, and roll angles.

[0011] Next, the sub-voxel positions and physical normal vectors of feature points are analyzed based on the network output results. For position analysis, firstly, response extrema are searched in the heatmap to determine integer coordinates. Then, a local window is defined centered on these coordinates. Using the predicted heatmap values ​​within the window as weights, a weighted average operation is performed on the voxel coordinates within the window to obtain sub-voxel coordinates. These sub-voxel coordinates are then mapped to physical space using the affine transformation matrix in the image metadata. For normal vector analysis, the pose parameters at the corresponding positions are read, a rotation matrix is ​​constructed according to a preset rotation order, and the preset reference axis vector is transformed into a local normal vector. To correct the influence of anisotropic scaling on the normal vector direction, the rotation and scaling components in the image affine transformation matrix are extracted, their inverse transpose is calculated, and this inverse transpose is used to transform and normalize the local normal vectors, ultimately obtaining the tangent plane normal vector in physical space.

[0012] Subsequently, feature point matching relationships are established. Since the edge feature points extracted by the network correspond to anatomical locations with fixed semantics on the surface of human organs or tissues, and each type of feature point is defined as unique in a single image, feature points with the same semantic category in the source data and target data directly form matching point pairs, without the need to perform similarity calculations based on image grayscale or local texture descriptors.

[0013] Finally, an objective function based on the distance from a point to the tangent plane is constructed, and the rigid body transformation parameters are solved. The objective function is constructed based on the following geometric constraint: after rigid body transformation, the edge feature points in the source data should lie on the tangent plane of the corresponding feature points in the target data. The specific calculation process is as follows: calculate the projection distance of the transformed source feature points on the normal direction of the tangent plane of the target feature points, and sum the squares of the projection distances of all matching point pairs as the optimization objective. To solve for the rotation matrix and translation vector corresponding to this objective function, a Gauss-Newton iterative algorithm or conjugate gradient method based on Lie algebra parameterization is adopted. This algorithm parameterizes the rotation matrix as a Lie algebra vector or uses the conjugate gradient method. In each iteration step, a weighted linear equation system is constructed, the Jacobian matrix and residual vector are calculated, the optimal perturbation amount is obtained by solving the linearized normal equation, and the rigid body transformation parameters are updated using exponential mapping operations. When the numerical amplitude of the perturbation amount is less than a preset threshold or the number of iterations reaches the upper limit, convergence is determined, and the final calibration result is output.

[0014] A second aspect of the present invention provides a computer device.

[0015] The device includes a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the method described in the first aspect above. The device is capable of automatically loading medical volume data, performing neural network inference, geometric feature analysis, iterative optimization, and image resampling tasks, ultimately generating registered volume data that is spatially strictly aligned with the target data.

[0016] The technical solution of this invention transforms the registration problem into a geometric alignment problem between points and tangent planes by directly predicting the normal vector of the tangent plane of anatomical landmarks. Compared with traditional point-to-point registration, the point-to-tangent plane model allows source feature points to slide tangentially within the target tangent plane, thereby significantly reducing the sensitivity of the registration results to errors in feature point extraction positions. Simultaneously, by directly regressing the normal vector parameters using a deep learning network, the computational burden and sensitivity to local noise associated with relying on fitting a plane to a neighborhood point cloud, as in traditional methods, are avoided. Furthermore, the introduction of a Lie algebra optimization framework ensures the orthogonal constraint properties of the rotation matrix during iterative updates, improving the stability and convergence efficiency of numerical computation.

[0017] This invention provides a method for edge feature point extraction and registration for rigid body registration of three-dimensional volume data. It has the following beneficial effects:

[0018] 1. This invention employs a point-to-tangent plane geometric constraint to construct the objective function. Compared to traditional point-to-point registration constraints, this method allows source feature points to slide tangentially within the tangent plane of the target feature points. This geometric mechanism reduces the dependence of the registration process on the accuracy of feature point extraction positions. Even if the feature points extracted by the neural network have a slight tangential offset on the anatomical surface, it will not produce a large distance penalty, thereby improving the fault tolerance of the registration algorithm for feature extraction errors.

[0019] 2. This invention utilizes independent branches of a fully convolutional neural network to directly predict the pose parameters of the tangent plane. Compared to traditional methods that search neighborhood point clouds and fit plane normal vectors, this direct prediction method avoids interference from local image noise in normal vector calculation, while also saving the computational overhead of complex neighborhood search and least squares fitting. End-to-end prediction using high-level semantic features extracted by the network ensures accurate geometric normal information is obtained even when image edges are blurred or artifacts are present.

[0020] 3. This invention introduces a Gauss-Newton iterative algorithm based on Lie algebra parameterization to solve for rigid body transformation parameters. By mapping the rotation matrix to the Lie algebra space for optimization, this method avoids the gimbaling problem that may arise from Euler angle representation and naturally maintains the orthogonality of the rotation matrix during iterative updates. This approach transforms complex nonlinear constraint optimization into an optimization problem on a linear manifold, effectively improving the stability of numerical computation and accelerating the convergence speed of registration iterations. Attached Figure Description

[0021] Figure 1 This is a flowchart of the edge feature point extraction and registration method for rigid body registration of 3D volume data;

[0022] Figure 2 A schematic diagram of the UNet network structure for extracting edge feature points;

[0023] Figure 3 A schematic diagram of edge feature points and tangent planes;

[0024] Figure 4 This is a geometric diagram showing the distance from a point to a tangent plane. Detailed Implementation

[0025] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0026] Please see the appendix Figure 1 - Appendix Figure 4 This invention provides a method for edge feature point extraction and registration for rigid body registration of three-dimensional volume data. It can be applied to a computer device containing a processor and a memory, wherein the processor is configured to execute a computer program stored in the memory to implement the following three-dimensional medical volume data processing logic. The method processes three-dimensional medical volume data stored in voxel matrix form, and the specific execution flow is as follows:

[0027] Step S1: Acquire the first and second 3D volume data to be registered. The first 3D volume data is defined as the source data, and the second 3D volume data is defined as the target data. Both the first and second 3D volume data are medical image data conforming to the DICOM standard, and their spatial positions are defined in the DICOM patient coordinate system. The DICOM patient coordinate system is based on the patient's anatomical structure, defining the X-axis pointing to the patient's left side, the Y-axis pointing to the patient's posterior side, and the Z-axis pointing to the patient's head side to ensure consistency of spatial geometric description across different scanning devices.

[0028] Step S2: Input the first 3D volume data and the second 3D volume data into the pre-trained UNet network respectively, and output the corresponding edge feature point heatmap and tangent plane pose parameter map.

[0029] The UNet network is constructed as a fully convolutional neural network, with its output containing two parallel and independent regression branches:

[0030] Sub-step S21: Generate an edge feature point heatmap through the first output branch. The edge feature point heatmap is a K-channel four-dimensional tensor, where K is the preset number of edge feature point categories. Each channel corresponds to a type of organ edge feature point with fixed anatomical semantics (e.g., lung apex, specific vertebral edge point), and this feature point is unique in anatomical structure. The voxel value in the heatmap represents the confidence probability that a feature point of the corresponding category exists at that location.

[0031] Sub-step S22: Generate the tangent plane attitude parameter map through the second output branch. The tangent plane attitude parameter map is a 3-channel four-dimensional tensor, and the values ​​of the three channels represent the pitch angle in the Eulerian system. ), yaw angle ) and roll angle (Roll, The attitude parameters are used to characterize the local surface geometry orientation at the corresponding voxel location.

[0032] Step S3: For each edge feature point category, analyze the geometric features and construct the tangent plane. For the k-th edge feature point category ( Perform the following processing:

[0033] Sub-step S31: Perform a full-image maximum search in the k-th channel of the edge feature point heatmap, and select the voxel coordinates with the largest response value as the unique edge feature point location of this category in three-dimensional space. .

[0034] Sub-step S32, coordinates of the attitude parameter diagram of the tangential plane At that location, read the corresponding pitch angle, yaw angle, and roll angle.

[0035] Sub-step S33: Construct the rotation matrix according to the preset rotation order of ZYX. Specifically, the rotation matrix Rotation about the Z-axis Rotation around the Y-axis Rotation around the X-axis The matrix is ​​obtained by multiplying the matrices.

[0036] Sub-step S34: Select a preset reference axis vector Using the rotation matrix The reference axis vector is transformed to obtain the normal vector of the tangent plane. The tangent plane is composed of the coordinates of the feature points. and its normal vector It is uniquely determined that its geometric equation satisfies .

[0037] Step S4: Establish a one-to-one matching relationship between the edge feature points of the first 3D volume data and the second 3D volume data based on the edge feature point categories. Since the edge feature point categories are defined based on anatomical semantics and are strictly unique, the k-th type edge feature point in the first 3D volume data directly forms a matching point pair with the k-th type edge feature point in the second 3D volume data, without the need to perform similarity calculation based on image grayscale or local descriptors.

[0038] Step S5: Construct an objective function based on the distance from the point to the tangent plane, and solve for the rigid body transformation parameters from the first three-dimensional volume data to the second three-dimensional volume data.

[0039] The rigid body transformation parameters include the rotation matrix. The translation vector is t. The objective function is configured as follows: calculate the distance from the edge feature points of the first 3D volume data, after rigid body transformation, projected onto the normal direction of the tangent plane of the corresponding edge feature points of the second 3D volume data, and sum the squares of this distance. The specific objective function is... Expressed as:

[0040]

[0041] in These are the coordinates of the k-th type of edge feature point in the first 3D volume data. It is the normal vector of its tangent plane; These are the coordinates of the k-th type of edge feature point in the second 3D volume data. It is the normal vector of its tangent plane; Let the rotation matrix be the rigid body transformation matrix from the first 3D volume data to the second 3D volume data. This is the rigid body transformation translation vector.

[0042] To ensure that the anatomical semantics of edge feature points and their tangent plane pose parameters are consistent across devices and scan sequences, this embodiment of the invention uniformly adopts the patient coordinate system defined by the DICOM standard for medical digital imaging and communication as the reference system for three-dimensional spatial calculation.

[0043] The DICOM patient coordinate system is a Cartesian coordinate system based on the anatomical orientation of the scanned object, i.e., the patient. In the DICOM patient coordinate system, the positive direction of the X-axis is defined as pointing to the left side of the patient; the positive direction of the Y-axis is defined as pointing to the back of the patient; and the positive direction of the Z-axis is defined as pointing to the head side of the patient. Therefore, the DICOM patient coordinate system is often referred to as the LPS coordinate system in the field of medical imaging. The unit of measurement for length in the DICOM patient coordinate system is defined as millimeters (mm).

[0044] The necessity of using the DICOM patient coordinate system lies in eliminating spatial differences introduced into medical image data during acquisition due to variations in scanning orientation or patient position. Medical image data is typically stored as a three-dimensional voxel matrix, which has its own image array coordinate system (indexed coordinate system). The axis of the image array coordinate system depends on the acquisition order and reconstruction direction of the scan slices. For example, transverse, coronal, or sagittal scans can cause the correspondence between the image array coordinate system and the anatomical orientation to be interchanged or flipped. If the positions of edge feature points or the orientation of the tangent plane are directly defined in the image array coordinate system, edge feature points with the same anatomical meaning will have completely different coordinate values ​​and normal vector directions in different scan data, thus compromising the uniqueness and registrationability of edge feature point categories.

[0045] This invention constructs an affine transformation matrix from the image array coordinate system to the DICOM patient coordinate system by reading the image position attribute, image orientation attribute, and pixel spacing attribute from the DICOM data file header. This affine transformation matrix incorporates the rotation direction cosine, the voxel physical size scaling factor, and the spatial origin translation vector. During UNet network inference and subsequent registration calculations, all spatial coordinate vectors involved refer to the physical coordinates in the DICOM patient coordinate system after transformation by the aforementioned affine transformation matrix.

[0046] Based on the DICOM patient coordinate system, this embodiment of the invention provides explicit physical definitions for the three angular components in the tangential plane attitude parameter diagram to support the subsequent construction of the rotation matrix. The angle values ​​output by the tangential plane attitude parameter diagram are expressed in radians, and the rotation direction follows the right-hand screw rule. The specific definitions are as follows:

[0047] The pitch angle output from the tangent plane attitude parameter diagram ( ) is defined as the angle of rotation about the X-axis (i.e., the left-right axis) of the DICOM patient coordinate system;

[0048] The yaw angle output from the tangential plane attitude parameter diagram ( ) is defined as the angle of rotation about the Y-axis (i.e., the front-to-back axis) of the DICOM patient coordinate system;

[0049] The roll angle output from the tangential plane attitude parameter graph ( ) is defined as the angle of rotation around the Z-axis (i.e., the head-to-foot axis) of the DICOM patient coordinate system.

[0050] The pose definition method based on the anatomical coordinate system ensures that the predicted cutting plane normal vector always maintains geometric semantic consistency with the patient's anatomical structure, regardless of the slice orientation of the original data, thus satisfying the prerequisite for establishing a one-to-one correspondence between feature points.

[0051] After establishing a unified DICOM patient coordinate system, this embodiment of the invention further clarifies the specific definition and selection criteria of the core geometric primitives used for registration calculation, namely the K edge feature point categories.

[0052] The edge feature point category refers to organ surface landmarks that have a clear name, fixed location, and can be stably identified across different individuals or scanning modalities in anatomical structures. Unlike conventional image feature points, the edge feature point category has clear medical semantic attributes. To ensure sufficient disclosure, taking the registration application of chest CT image data as an example, the edge feature point category may specifically include, but is not limited to: the highest anatomical point of the left upper lobe of the lung, the highest anatomical point of the right upper lobe of the lung, the center point of the tracheal carina, the center point of the sternal notch, and the tangent point at the top of the aortic arch. For the registration application of orthopedic image data, the edge feature point category may specifically include: the geometric center point of the femoral head, the lateral point of the greater trochanter of the femur, the lower edge of the ischial tuberosity, etc. Each of the aforementioned edge feature point categories corresponds to a fixed physiological structural position in the human anatomical atlas.

[0053] The edge feature point category is defined as strictly unique in each 3D volume data. Strict uniqueness means that for any given edge feature point category, there exists only one voxel coordinate that forms a mapping relationship with it in the voxel set of a single 3D volume data, without ambiguity or multiple solutions. To ensure this strict uniqueness in algorithm implementation, embodiments of the present invention use geometric extremum constraints or anatomical topological constraints to define feature points. For example, the highest anatomical point of the left upper lobe of the lung is defined as the voxel with the largest Z-axis coordinate value within the left lung mask region; the outermost point of the greater trochanter of the femur is defined as the voxel with the smallest coordinate value in the negative X-axis direction of the femur region. This definition method based on geometric extremum ensures that regardless of the resolution of the 3D volume data, the edge feature point category always mathematically converges to a unique spatial location.

[0054] Furthermore, the edge feature point category is specifically defined as an edge feature point, meaning that the physical location corresponding to the edge feature point category must be located on the surface boundary of a human organ or tissue, rather than a filled area inside the organ. The technical significance of this limitation is that points located on the surface of an object naturally possess local surface geometric properties, i.e., they have definite tangent planes and normal vectors. If the feature point is located inside the tissue, such as the center of the liver, a physically meaningful tangent plane cannot be defined, thus making it impossible to apply the registration algorithm based on the distance from the point to the tangent plane described in this embodiment of the invention. Therefore, each edge feature point category, in the corresponding three-dimensional volume data, not only has unique three-dimensional coordinates but is also associated with a unique tangent plane that is tangent to the surface of the organ at that location.

[0055] The preset number of edge feature point categories, K, is typically greater than or equal to 3, and the selected edge feature point categories should be non-collinearly distributed in three-dimensional space as much as possible. By selecting edge feature point categories distributed in different orientations of the anatomical structure, such as top, bottom, left, right, front, and back, the six degrees of freedom of rigid body transformation can be constrained to the greatest extent, namely three translational degrees of freedom and three rotational degrees of freedom, avoiding the problem of unstable rotation angle calculation caused by overly concentrated feature point distribution. During the training phase, the definitions of the K edge feature point categories are fixed in the manually labeled rules, and their anatomical semantic knowledge is mapped to the weight parameters of the UNet network through supervised learning, enabling the UNet network to automatically and accurately locate anatomical feature points that conform to the above definitions during the inference phase.

[0056] In step S2, the edge feature point extraction network first performs deep feature extraction on the input 3D volume data through a shared backbone network to generate a common feature basis, and then branches to the first output branch in step S21 and the second output branch in step S22.

[0057] The edge feature point extraction network adopts a fully convolutional 3D neural network architecture, and the 3D volume data received at its input is formatted into a dimension of _____. The tensor, where 1 represents a single-channel grayscale input. These represent the depth, height, and width of the 3D volume data, respectively. Before executing steps S21 and S22, the edge feature point extraction network sequentially performs the following backbone feature extraction sub-steps:

[0058] Sub-step S201: Perform multi-level feature encoding.

[0059] The feature encoding path is configured to progressively downsample and abstract the input 3D volume data. This path comprises four cascaded encoding stages, each sequentially performing convolutional feature extraction and downsampling operations:

[0060] 1. Convolutional Feature Extraction: Each encoding stage contains two consecutive convolutional blocks. Each convolutional block consists of a feature vector of size [missing information]. The system consists of a 3D convolutional layer, a normalization layer, and a non-linear activation function layer connected in series. The stride and padding of the 3D convolutional layer are set to 1 to maintain the feature map size. The normalization layer uses instance normalization to adapt to training scenarios with smaller medical images. The non-linear activation function layer uses a modified linear unit with leakage, and its negative half-axis slope is set to 0.01 to prevent gradient vanishing.

[0061] 2. Downsampling operation: After two convolutional blocks, a 3D max pooling operation with a stride of 2 is performed. This operation downsamples the feature map... The dimensions in each dimension are halved, while the number of feature channels is doubled. Through four encoding stages, the number of channels in the feature map increases sequentially from the initial 64 to 128, 256, and 512, forming a multi-scale encoded feature pyramid.

[0062] Sub-step S202: Perform bottleneck layer feature abstraction.

[0063] The bottleneck layer is connected to the end of the feature encoding path and is used to process the feature map with the lowest spatial resolution but the highest semantic dimension. The bottleneck layer contains two consecutive 3×3×3 convolutional blocks, with the same structure as above, and its feature channel number is set to 1024. The bottleneck layer does not perform downsampling operations; its function is to integrate global contextual information and establish a high-level anatomical semantic representation.

[0064] Sub-step S203: Perform feature decoding and skip connection fusion.

[0065] The feature decoding path, located after the bottleneck layer, comprises four cascaded decoding stages, aiming to progressively restore the spatial resolution of the feature map and fuse shallow detail information:

[0066] 1. Upsampling operation: In each decoding stage, the spatial size of the input feature map is first increased by 2 times through 3D transposed convolution or trilinear interpolation, and the number of channels is halved.

[0067] 2. Skip Connection Fusion: The upsampled feature map is concatenated with the shallow feature map generated in the encoding stage corresponding to the spatial resolution in sub-step S201. The shallow feature map retains high-frequency spatial information such as edges and textures of the original image, and the concatenation operation compensates for the spatial localization accuracy lost by the deep network.

[0068] 3. Convolutional Feature Reassembly: Two consecutive 3×3×3 convolutional blocks are applied to the concatenated feature map to fuse deep semantic features with shallow spatial features.

[0069] After processing in sub-step S203, the edge feature point extraction network outputs a common feature map restored to the original input resolution (D, H, W). This common feature map is then simultaneously fed into the aforementioned first output branch (corresponding to step S21) and the second output branch (corresponding to step S22), respectively, to map and generate the edge feature point heatmap and the tangent plane pose parameter map.

[0070] For the first output branch in step S2 that generates the heatmap of edge feature points, corresponding to step S21 above, this embodiment of the invention designs a dedicated regression prediction structure and defines a supervision signal generation rule based on a three-dimensional Gaussian distribution to ensure that the network can accurately learn the spatial distribution of sparse feature points. The specific construction and processing logic of the first output branch is as follows:

[0071] Step S211: Construct the heat map prediction layer structure.

[0072] The first output branch directly receives the common feature map output by the backbone of the edge feature point extraction network. To map high-dimensional semantic features to a probability distribution for a specific category, the first output branch is initially configured with a... The three-dimensional convolutional layer. 3D convolutional layers contain There are n convolutional kernels, among which This corresponds to the preset number of edge feature point categories. Through this convolution operation, the channel dimension of the common feature map is compressed and mapped from the feature dimension to the category dimension. (channels), while maintaining the same three-dimensional spatial resolution.

[0073] In the Following the 3D convolutional layer, a cascaded Sigmoid non-linear activation function layer is applied. This Sigmoid non-linear activation function layer maps the real-number domain values ​​of the convolutional output to the (0,1) interval, thus giving the output tensor a physical meaning of probability confidence. After this processing, the first output branch finally outputs a tensor with dimension 1. The edge feature point heatmap. In the edge feature point heatmap, the first... Each voxel value in the channel represents the voxel's position at that spatial coordinate, indicating that it belongs to the th channel. Confidence probability of edge feature points.

[0074] Step S212: Define the rules for generating heatmap monitoring signals.

[0075] To train the first output branch, a realistic label heatmap with the same size as the model output needs to be constructed. This embodiment of the invention does not use binary one-hot encoding as labels, but instead uses a three-dimensional Gaussian kernel function to generate smooth soft labels, thereby alleviating the spatial uncertainty of feature point definition in medical images and providing a denser gradient backpropagation signal. For the first... For edge feature points, let their manually labeled real physical coordinates be... In the first section of the real label heatmap In each channel, arbitrary voxel coordinates The value at the location The following three-dimensional Gaussian distribution formula is used for calculation:

[0076]

[0077] The symbols are defined as follows:

[0078] Represents the Euclidean norm, used to calculate voxel coordinates. Coordinates of the true feature points The spatial Euclidean distance between them; The standard deviation of the Gaussian kernel is used to control the dispersion of the response region in the heatmap. In a specific embodiment, the... The value is set to a preset voxel distance to cover the local neighborhood around the feature point; Represented by natural constant An exponential function with base 0.

[0079] Step S213: Configure the multi-task joint loss function. To achieve end-to-end collaborative training of edge feature point positions and tangent plane poses, this embodiment of the invention constructs a joint loss function that includes a heatmap regression error term and a pose parameter error term. During network training, the joint loss function is minimized... To synchronize and update network weights. The joint loss function... The specific calculation formula is as follows:

[0080]

[0081] The physical meaning and calculation rules of each symbol in the formula are defined as follows:

[0082] : Represents the total loss value used for backpropagation; : Represents the total number of preset edge feature point categories; : Represents the voxel space coordinate index in 3D volume data, summation symbol This indicates that all voxel positions are traversed; : indicates the first output branch prediction obtained by the first branch. Class edge feature points in coordinates The confidence level of the heatmap at the location; : indicates the generation rule obtained according to step S212. The true label heatmap values ​​of class edge feature points; : This represents the voxel-level loss function that measures the difference in heatmap predictions, specifically using the mean squared error function or the focal loss function; : indicates that the network is in the . The predicted tangent plane pose parameter vector is obtained by sampling the true coordinates of the edge feature points; the sampling operation is performed on the predicted pose parameter map by sampling the true coordinates. Obtained by index extraction; : indicates the first Manually labeled true tangent plane pose parameter vectors corresponding to class edge feature points; : Represents the vector distance metric norm, specifically using the L1 norm (absolute error) or the L2 norm; and : These represent the balancing weight coefficients for the heatmap regression loss term and the attitude parameter loss term, respectively, used to balance the contribution ratio of the two tasks in the gradient descent process.

[0083] Through the design of the joint loss function described above, the edge feature point extraction network is not only constrained to produce high response values ​​at the correct locations (by the first term) (Control), while being forced to predict the correct geometric tangent plane pose at the feature point location (by the second term) (Control), thereby achieving decoupled regression of position and orientation.

[0084] Parallel to the first output branch described in step S21, this embodiment of the invention configures a second output branch at the end of the edge feature point extraction network, corresponding to the tangent plane pose prediction branch in step S22. Specifically, the second output branch is configured to decode the local geometric orientation information of the edge feature points from the common feature map to support subsequent rigid body registration calculations based on the tangent plane. The specific construction and processing logic of the second output branch is as follows:

[0085] Step S221: Construct the attitude parameter prediction layer structure.

[0086] The second output branch receives the common feature map output by the backbone of the edge feature point extraction network. To map the high-dimensional semantic features to rotational pose parameters in three-dimensional space, the second output branch is configured with a... The three-dimensional convolutional layer. The 3D convolutional layer contains three convolutional kernels, the number of which strictly corresponds to the three Euler angle degrees of freedom required to describe the rotation of a 3D rigid body. This convolutional operation compresses and maps the channel dimension of the common feature map to 3 channels while maintaining the 3D spatial resolution.

[0087] In the Following the 3D convolutional layer, a hyperbolic tangent activation function layer is cascaded. This hyperbolic tangent activation function layer forces the output value of the convolution to be within the range of [-1, 1]. To cover the complete rotation angle range in 3D space, this embodiment further includes a scaling operation, multiplying the output value of the hyperbolic tangent activation function layer by the constant pi. After the above processing, the second output branch finally outputs a dimension of... The attitude parameter diagram of the tangent plane.

[0088] Step S222: Define the channel physical semantics of the attitude parameter map.

[0089] In the tangential plane attitude parameter diagram, the values ​​of the three channels each have a clear geometric and physical definition to maintain consistency with the aforementioned definition of the DICOM patient coordinate system. In this embodiment of the invention, the initial tangential plane normal vector is defined as a unit vector pointing in the positive Z-axis direction of the DICOM patient coordinate system. The three angular components output from the tangent plane attitude parameter map are used to construct a rotation matrix, rotating the initial tangent plane normal vector to the target tangent plane normal vector. The specific definitions are as follows:

[0090] The value output by the first channel of the tangential plane attitude parameter map corresponds to the pitch angle. The value represents the radian value of rotation about the X-axis of the DICOM patient coordinate system, and its range is:

[0091] The value output by the second channel of the tangential plane attitude parameter graph corresponds to the yaw angle. The value represents the radian value of rotation about the Y-axis of the DICOM patient coordinate system, and its range is:

[0092] The value output by the third channel of the tangential plane attitude parameter graph corresponds to the roll angle (Roll, The value represents the radian value of rotation about the Z-axis of the DICOM patient coordinate system, and its range is: .

[0093] To eliminate the non-commutative ambiguity of rotation operations, embodiments of the present invention stipulate that the rotation sequence follows a fixed order of either intrinsic or extrinsic rotations, for example, following... The matrix multiplication order. Through this dense prediction method, the edge feature point extraction network outputs a set of pose parameters for each voxel position in the input volume data to determine the local tangent plane normal vector.

[0094] Step S223: Perform sparse sampling supervised training.

[0095] Unlike the first output branch, which employs dense supervision across the entire image, the second output branch is trained using a sparse supervision strategy. This is because the surface tangent planes of the anatomical structure only have a clear physical definition at the actual edge feature point locations.

[0096] During training, manually labeled data was used. The true physical coordinates of edge feature points The voxel index coordinates are obtained by rounding, and these coordinates are used as indices to extract the predicted attitude vector at the corresponding position from the tangent plane attitude parameter map output by the second output branch. The predicted attitude vector It includes the pitch angle, yaw angle, and roll angle predicted for that coordinate position.

[0097] Subsequently, the predicted pose vector Compared with manually labeled true pose vectors Substitute the values ​​into the joint loss function formula described in step S213 above to calculate the error. During backpropagation, gradient information is obtained only through real physical coordinates. The location of the voxel is fed back to the edge feature point extraction network, allowing the network to update weights only for the local geometric surface orientation at specific anatomical landmarks, while ignoring pose predictions for non-feature point regions. This sparse sampling mechanism effectively reduces the learning difficulty and decouples the position regression task from the pose estimation task.

[0098] To construct a supervisory signal for driving the edge feature point extraction network to optimize parameters, this embodiment of the invention formulates a detailed Gaussian heatmap label generation strategy. This strategy aims to transform discrete anatomical landmark coordinates into a continuous, differentiable probability distribution field, overcoming subtle errors inherent in manual labeling and providing spatial gradient guidance. Since the three-dimensional Gaussian distribution calculation relationship of the heatmap values ​​has been clearly defined in step S212, this section focuses on elaborating the steps for transforming this mathematical definition into a computer-executable engineering implementation, specifically including the following process:

[0099] Step S301: Perform voxel space mapping of edge feature point coordinates.

[0100] Raw annotation information for medical image data is typically stored in physical coordinates within the DICOM patient coordinate system. To generate labels aligned with the network's output feature map, these physical coordinates must first be mapped to voxel index coordinates in the image array space.

[0101] This invention reads data recorded in medical image metadata. The affine transformation matrix M describes the transformation relationship from the voxel index space to the physical space. Therefore, for the ... Physical coordinates of class edge feature points In this embodiment of the invention, the inverse matrix of the affine transformation matrix is ​​calculated. This is used to achieve the inverse mapping. The calculation process employs a homogeneous coordinate transformation, that is, first constructing an expanded coordinate vector. Then multiply by the inverse matrix on the left. voxel index coordinates are obtained It should be noted that the calculated voxel index coordinates Preserve floating-point precision and do not perform rounding operations to retain sub-voxel level position information to the greatest extent possible.

[0102] Step S302: Implement coordinate adaptation for the target resolution.

[0103] Considering that edge feature point extraction networks typically involve downsampling operations, the spatial resolution of the network output layer is often lower than the resolution of the original input image. To ensure the consistency between the supervision signal and the network output in the spatial dimension, this embodiment of the invention obtains the bus sampling rate of the edge feature point extraction network. (For example .

[0104] The voxel index coordinates obtained in step S31 Divide by the downsampling factor To obtain target coordinates adapted to the network output layer resolution. Subsequent heatmap generation operations are all based on this target coordinate. Performed in the scaled grid space.

[0105] Step S303: Perform Gaussian kernel filling based on local truncation.

[0106] To improve computational efficiency and reduce memory usage, this embodiment of the invention employs a local truncation strategy instead of full-image computation. For each edge feature point category, this embodiment of the invention only performs local truncation calculations using the target coordinates. Calculate the heatmap values ​​within the local cubic region centered on the heatmap.

[0107] The side length range of the local cubic region is set to... The cutoff radius Set as Gaussian kernel standard deviation The value is three times that of the local cube region. Within the local cube region, the value of each voxel is calculated strictly according to the three-dimensional Gaussian distribution formula defined in step S212 above; for voxels located outside the local cube region, their heatmap values ​​are directly assigned to 0. This processing method utilizes the rapid decay characteristic of the Gaussian function, significantly reducing computational complexity while ensuring numerical accuracy.

[0108] Step S304: Handling missing feature points and tensor assembly.

[0109] In actual clinical data, due to limitations in scanning range, some anatomical landmarks may not be included in the image data. To address this issue, this embodiment of the invention assigns a value of 0 to the entire heatmap channel corresponding to that category and generates a corresponding mask tensor. This mask tensor is used to shield the gradient backpropagation of that channel when calculating the loss function, preventing the network from learning incorrect background features.

[0110] After completing all After generating heatmaps for each category, this embodiment of the invention stacks the individual channel heatmaps along the channel dimension to assemble them into a heatmap with dimension [missing information]. The complete real-world label heatmap tensor is used for subsequent network training.

[0111] To facilitate the training of the tangential plane pose prediction branch constructed in step S22, this embodiment of the invention constructs a tangential plane pose label generation process based on the geometric characteristics of anatomical structure surfaces, and designs a sparse supervision mechanism to improve training convergence efficiency. The specific implementation of the tangential plane pose label and sparse supervision includes the following steps:

[0112] Step S305: Extract the true surface normal vector at the anatomical landmark.

[0113] To obtain the true label of the tangent plane pose, it is necessary to acquire the geometric normal vectors of the edge feature points at the anatomical structure surface. This embodiment of the invention utilizes an anatomical structure segmentation mask or a surface triangular mesh model corresponding to the medical image data for calculation. When using an anatomical structure segmentation mask, this embodiment performs spatial gradient calculation on the binarized anatomical structure segmentation mask. For the physical location of the kk-th class of edge feature points, the direction of the spatial gradient vector of the anatomical structure segmentation mask at that location indicates the direction of the surface normal at that location. When using a surface triangular mesh model, triangular patches containing the kk-th class of edge feature points are located, and the geometric normal vectors of these triangular patches are used as the true surface normal vectors of the edge feature points. The true surface normal vectors are then normalized to a magnitude of 1.

[0114] Step S306: Convert the real surface normal vector to Euler angle labels.

[0115] Since the second output branch in step S22 is configured to regress and predict Euler angle parameters, it is necessary to convert the true surface normal vector into the corresponding Euler angle form. Based on the geometric definition of the three-dimensional Cartesian coordinate system, the component values ​​of the true surface normal vector on each coordinate axis directly reflect the inclination of the tangent plane relative to the coordinate plane. In this embodiment of the invention, based on the projection components of the true surface normal vector on the X, Y, and Z axes of the DICOM patient coordinate system, the pitch and yaw angle values ​​required to rotate the initial reference vector to coincide with the true surface normal vector are calculated analytically through reverse geometric projection. For the roll angle degree of freedom of rotation around the normal vector, since the rotation state within the tangent plane cannot be uniquely determined solely by the normal vector, this embodiment of the invention introduces a zero roll constraint. The zero roll constraint fixes the true label value of the roll angle to 0, thereby establishing a unique mapping relationship from the true surface normal vector to the Euler angle attitude parameters.

[0116] Step S307: Construct a sparse supervised mask and gradient blocking.

[0117] Considering that the tangent plane pose has clear physical meaning only at feature points on the organ surface, and is meaningless in the background region or inside the organ, this embodiment of the invention implements a sparse supervision strategy during training. The specific implementation of the sparse supervision strategy is as follows: a binary weight mask tensor with the same size as the network output is constructed. In the binary weight mask tensor, the voxel weight is set to 1 only at the voxel coordinates corresponding to the true positions of edge feature points; and the voxel weight is set to 0 at all other voxel coordinates. When calculating the pose parameter loss term described in step S213, the pose parameter map predicted by the network is multiplied element-wise with the binary weight mask tensor. This operation is logically equivalent to extracting only the prediction error of the feature point positions for backpropagation. Through this sparse supervision mechanism, during the network backpropagation stage, the weight gradient of the pose prediction branch is driven only by the error of the feature point positions. The random noise prediction values ​​generated in the background region do not produce gradient signals, thereby achieving gradient blocking of invalid regions and significantly improving the convergence speed of the network on sparse feature learning tasks.

[0118] Step S308: Perform geometrically consistent data augmentation transformation.

[0119] During network training, to prevent overfitting, spatial data augmentation operations such as random rotation and flipping are typically performed on the input 3D volume data. This embodiment of the invention emphasizes the geometric consistency between pose labels and spatial transformations. When a random rotation matrix is ​​applied to the input 3D volume data, the ground truth surface normal vector must be simultaneously left-multiplied by the random rotation matrix to obtain a new transformed normal vector. Subsequently, the transformation process described in step S36 is repeated, and the corresponding transformed pose parameter labels are recalculated using the new normal vector. This dynamic label update mechanism ensures that regardless of how the input image is rotated, the pose parameters learned by the network always maintain the correct geometric orientation relative to the current transformed image coordinate system, achieving equivariance of feature extraction with spatial geometric transformations.

[0120] To achieve simultaneous optimization of the location regression and pose estimation tasks in the edge feature point extraction network, this embodiment of the invention constructs a multi-task joint loss function. This multi-task joint loss function measures the difference between the ground truth labels generated in the preceding steps and the network prediction results, and drives network weight updates through a backpropagation algorithm. The specific implementation of the loss function construction includes the following steps:

[0121] Step S309: Construct a joint objective function for multiple tasks.

[0122] The edge feature point extraction network is a multi-branch output architecture. To balance the accuracy of heatmap regression in the first output branch with the accuracy of pose parameter estimation in the second output branch, this embodiment defines the total loss function as a weighted sum of the position regression loss term and the pose estimation loss term. The formula for calculating the total loss function is defined as follows:

[0123]

[0124] The symbols are defined as follows: This represents the location regression loss term used to measure the accuracy of heatmap predictions; This represents the attitude estimation loss term used to measure the accuracy of tangential plane attitude prediction; The weighting coefficients of the location regression loss term are used to normalize the gradient magnitude for the heatmap regression task. The weighting coefficients of the attitude estimation loss term are used to adjust the contribution of the attitude estimation task to the overall optimization process. In a specific embodiment, the weighting coefficients of the position regression loss term... With the weighting coefficients of the attitude estimation loss term The numerical settings are balanced based on the order of magnitude of the losses of the two tasks.

[0125] Step S310: Calculate the heatmap location regression loss term.

[0126] The location regression loss term is configured to measure, at the voxel level, the difference between the predicted heatmap output by the first output branch and the true Gaussian heatmap generated in steps S31 to S34. This embodiment of the invention uses mean squared error as the basic metric function. The specific formula for calculating the location regression loss term is as follows:

[0127]

[0128] The symbols are defined as follows: This represents the total number of categories of edge feature points; This represents the total number of voxels in the feature map; Indicates the first output branch prediction of the first branch. The first channel Confidence level of individual element location; Represents the first Gaussian heat map. The first channel The label value of the individual element location; the summation symbol This represents the summation of errors across all categories of channels and all spatial voxel positions. By minimizing this location regression loss term, the heatmap distribution output by the edge feature point extraction network will approximate the true Gaussian distribution, thereby achieving feature point localization.

[0129] Step S311: Calculate the attitude estimation loss term under sparse mask constraints.

[0130] The pose estimation loss term is configured to calculate the pose prediction error only for actual anatomical landmark locations. To achieve this, this embodiment of the invention utilizes the binary weighted mask tensor constructed in step S37 above to filter out predicted values ​​for invalid background regions. This embodiment of the invention uses the absolute value error of the L1 norm as the basic metric for pose regression. The specific formula for calculating the pose estimation loss term is as follows:

[0131]

[0132] The symbols are defined as follows: The binary weight mask tensor represents the first... The mask value corresponding to the position of the nth feature point, when the nth feature point is... Individual element corresponds to the first The value is 1 when the actual location of the feature point is known, and 0 otherwise. This represents the total number of valid feature points across all categories, i.e., the number of non-zero elements in the mask tensor. It represents a very small positive constant, used to prevent numerical calculation errors when the denominator is zero; , , These represent the network at the 1st and 2nd epochs, respectively. The predicted pitch, yaw, and roll angles for each position; , , These represent the actual Euler angle label values ​​calculated in the preceding steps; This represents absolute value operations.

[0133] It should be noted that, to avoid gradient calculation errors caused by angular periodicity, this embodiment of the invention normalizes all angle values ​​to a linear continuous interval when generating labels, ensuring that the L1 loss function can correctly measure angle differences. This is achieved by introducing a mask term. The above formula ensures that the angle prediction error only contributes to the network gradient when the feature point actually exists at its coordinate position, thus achieving sparse supervision.

[0134] Step S312: Perform backpropagation and optimizer iteration.

[0135] After calculating the total loss function, this embodiment of the invention uses an automatic differentiation mechanism to calculate the gradient of the total loss function with respect to all trainable parameters in the edge feature point extraction network. Subsequently, a stochastic gradient descent or adaptive moment estimation optimizer is used to update the network parameters based on the calculated gradient direction. During training, this embodiment of the invention also configures a learning rate decay strategy, gradually reducing the learning rate value as the number of training iterations increases to ensure that the network model can converge to the optimal solution in the later stages of training.

[0136] After the edge feature point extraction network is trained, this embodiment of the invention enters the inference phase. The main task of the inference phase is to use the trained network model to accurately extract the physical coordinates of edge feature points from the input medical image data. For the specific task of feature point coordinate extraction, this embodiment of the invention designs a cascaded processing flow from coarse localization to sub-voxel fine correction, specifically including the following steps:

[0137] Step S41: Obtain the network prediction heatmap and integer peak index.

[0138] The target medical image data to be processed is input into the trained edge feature point extraction network. The edge feature point extraction network performs a forward propagation operation and outputs a data set containing edge feature point data through the first output branch. The predicted heatmap tensor for each channel. The dimension of the predicted heatmap tensor is... ,in These correspond to the depth, height, and width dimensions of the feature map, respectively.

[0139] For the The feature point category, the computer processing unit in the corresponding _th ... A full scan is performed in each heatmap channel to find the voxel location with the largest response value. The location of the largest response value is determined as an integer-order peak coordinate. The integer-level peak coordinates represent the approximate positions of edge feature points on the discrete feature map grid.

[0140] To eliminate interference from image noise or missing feature points, this embodiment of the invention sets a confidence threshold. (For example If the response value at the integer-level peak coordinate is less than the confidence threshold... Then determine the first If a feature point of class 1 is not visible in the current image, stop processing that feature point. Subsequent calculations of class feature points; if the response value is greater than or equal to the confidence threshold. If so, continue with the subsequent refined positioning steps.

[0141] Step S42: Subvoxel coordinate correction based on the local centroid method.

[0142] Since the integer-level peak coordinates are limited by the discrete resolution of the feature map and cannot meet the requirements of high-precision surgical navigation, this invention uses the local centroid method to calculate sub-voxel-level coordinates. Define a local cube window centered on the cube. (For example, the size is) In the local cube window Within the system, the weighted average of spatial coordinates is calculated using the values ​​from the predicted heatmap as weights, thus obtaining sub-voxel level coordinates.

[0143] The subvoxel level coordinates The calculation formula is as follows:

[0144]

[0145] The symbols are defined as follows: Represents a local cube window An integer index coordinate vector of any voxel within the vector; Indicates the first Each heatmap channel is located at coordinates The predicted response value at the location; summation symbol This indicates that all voxels within the window are iterated and accumulated. Using the local centroid method, the distribution information within the peak neighborhood is used to interpolate the position, enabling the final extracted coordinates to overcome the limitations of grid resolution and achieve sub-voxel level positioning accuracy.

[0146] Step S43: Physical space coordinate mapping and restoration.

[0147] The output feature map of the edge feature point extraction network is downsampled and scaled compared to the original input image, and the currently obtained sub-voxel level coordinates... The image is still in the image array coordinate system. To obtain the true physical coordinates for surgical navigation, a spatial mapping transformation needs to be performed. This involves acquiring the raw spatial metadata of the medical image data, including the downsampling rate from the network preprocessing stage. and the affine transformation matrix corresponding to the image The affine transformation matrix The transformation relationship from image voxel space to DICOM patient physical coordinate system is described.

[0148] Using the downsampling factor and the affine transformation matrix The subvoxel level coordinates Converted to final physical feature point coordinates The calculation process first multiplies the sub-voxel level coordinates by the scaling factor. To restore the image to its original resolution scale, a homogeneous coordinate vector is then constructed and left-multiplied. affine transformation matrix .

[0149] The coordinates of the physical feature points The specific calculation formula is as follows:

[0150]

[0151] in, Indicates after scale restoration Homogeneous coordinate vector. Calculation result. The first three components constitute the first The precise three-dimensional spatial location of the edge feature points on the patient's anatomical structure will serve as the key input point for the subsequent rigid body registration algorithm.

[0152] After obtaining the physical coordinates of the edge feature points, in order to fully construct the joint position and pose constraints for surgical navigation, this embodiment of the invention continues to execute the parsing process of the tangential plane pose. This process uses the prediction results of the second output branch of the edge feature point extraction network to reconstruct the normal vector of the anatomical tangential plane passing through the edge feature points in physical space. The specific implementation of the tangential plane pose parsing includes the following steps:

[0153] Step S44: Extract the pose prediction parameters at the feature point locations.

[0154] The second output branch of the edge feature point extraction network outputs a pose parameter tensor. The spatial dimensions of the pose parameter tensor are the same as those of the heatmap tensor of the first output branch.

[0155] Using the integer-level peak coordinates of the k-th type feature point determined in step S41 above. As a spatial index, three values ​​are extracted from the corresponding coordinate positions of the attitude parameter tensor. These three values ​​correspond to the predicted pitch angle at the k-th type feature point. Predicted yaw angle And predicting the roll angle .

[0156] In practice, if the network output consists of multiple shared channels, it is read directly; if it consists of independent channels for each category, the channel group corresponding to the k-th category is read. This keypoint index-based reading method avoids the calculation of invalid pose information in the background region and effectively filters out noise interference from non-anatomical regions.

[0157] Step S45: Reconstruct the unit normal vector in the local coordinate system.

[0158] Obtained Euler angle parameters , , This describes the rotation relationship of the tangent plane normal vector relative to the image voxel coordinate system. To recover the vector form of the normal vector, this embodiment performs a reverse derivation based on the geometric mapping relationship defined in step S36 of the training phase. Based on the aforementioned zero roll angle constraint and the principle of spherical coordinate system transformation, this embodiment maps the predicted Euler angles back to the three-dimensional Cartesian coordinate system, constructing a unit normal vector in the image voxel space. The unit normal vector The specific calculation formula is as follows:

[0159]

[0160]

[0161]

[0162] The above formula is directly derived from the positive mapping relationship defined during the training phase, ensuring the consistency between the training objective and the inference result in terms of geometric definition. The calculated unit normal vector... It characterizes the orientation of the normal of the tangent plane in the image array space of the medical image.

[0163] Step S46: Perform physical space transformation of the normal vector.

[0164] Because medical imaging data typically exhibits anisotropic physical resolution and rotational deviation relative to the patient's coordinate system, the unit normal vector... This only represents the orientation in the image grid space. To obtain the true physical normal vectors guiding surgical instrument manipulation, they must be transformed to the DICOM patient physical coordinate system. This embodiment of the invention reads data from medical image metadata... Affine transformation matrix Since the normal vector is a direction vector, its transformation is not affected by the translation components of the coordinate system, but only by the rotation and scaling components. Therefore, the top left corner of the affine transformation matrix M is extracted. Submatrix, denoted as rotation and scaling matrix According to the principles of differential geometry, surface normal vectors are covariant vectors, and their transformation rules follow the inverse inverse form of the transformation matrix. This embodiment of the invention calculates the rotation-scaling matrix. inverse matrix And use the inverse matrix to apply the unit normal vector The transformation is performed, followed by a normalization operation, to obtain the final physical tangent plane normal vector. The physical tangent plane normal vector The calculation formula is as follows:

[0165]

[0166] in, This represents the L2 norm operation, used to ensure that the output vector is a unit vector. Up to this point, this embodiment of the invention has fully obtained the [number of] [vectors / properties]. Physical coordinates of edge feature points The result is obtained from step S43, along with the physical normal vector of the tangent plane passing through that point. The physical location coordinates They jointly defined a definite spatial plane equation This plane equation accurately characterizes the surgical cutting plane of the target anatomical structure, providing complete geometric constraints for subsequent surgical navigation registration and instrument positioning.

[0167] After successfully parsing the feature point set with tangent plane pose information from the medical image data in the aforementioned steps, this embodiment of the invention enters the rigid body registration solution stage. The core task of this stage is to calculate the rigid body transformation matrix between the patient's actual physical space coordinate system and the medical image coordinate system. This embodiment of the invention proposes a high-precision registration method based on a point-to-tangent plane distance metric. First, a feature matching relationship with explicit semantic constraints is established. The specific implementation includes the following steps:

[0168] Step S51: Construct the data structure for the source point set and the target tangent plane set.

[0169] This invention relates to geometric datasets in two independent coordinate systems. The first dataset is the source point set. Defined within the tracking coordinate system of the surgical navigation system. The computer processing unit receives discrete coordinate point data acquired by the user using a spatial tracking device on the actual anatomical structure surface of the patient. For the predefined... Each anatomical feature point category is used to obtain the included... A set of three-dimensional coordinate points .in, Indicates the first The actual measured coordinates of the feature points in the tracking coordinate system.

[0170] The second dataset is the target tangent plane set. Defined in the physical coordinate system of medical image data, it originates from the inference results of the aforementioned steps S43 and 546. For the same The edge feature point extraction network outputs the corresponding position coordinates in the medical image coordinate system for each anatomical feature point category. and tangent plane normal vector Therefore, the target tangent plane set is constructed as follows:

[0171]

[0172] Step S52: Establish feature point matching relationships based on semantic consistency.

[0173] To avoid the iterative search computation process required in traditional registration algorithms due to the lack of prior correspondence, this embodiment of the invention utilizes the semantic classification characteristics of the edge feature point extraction network to establish a deterministic index mapping mechanism.

[0174] Due to the source point set Each sampling point in All are associated with specific anatomical landmark category labels, and the target sectional plane set Each element in All are generated by the network in the corresponding number of The outputs in each feature channel have a unique correspondence in terms of anatomical semantics.

[0175] According to the category index of the present invention Directly build Matching pairs This semantically consistent forced matching mechanism ensures that each physical space sampling point precisely corresponds to the tangent plane model of the same anatomical location in the medical image space. For potentially invalid points, i.e., points whose confidence level is lower than the threshold in the aforementioned step S41, resulting in no output results, this embodiment of the invention performs an intersection filtering operation when constructing matching pairs, retaining only valid feature point indices that exist in both the source point set and the target tangent plane set for subsequent calculations.

[0176] After establishing the semantic matching relationship between the source point set and the target tangent plane set, this embodiment of the invention continues to perform the core step of rigid body registration, namely, constructing an objective function to describe the registration error. This embodiment uses a distance metric from a point to the tangent plane to construct the objective function. This metric reduces the impact of source point position deviations on registration convergence by ignoring the tangential displacement components within the tangent plane. The specific implementation of the objective function construction includes the following steps:

[0177] Step S53: Define rigid body transformation parameter variables.

[0178] To map the source point set in the surgical navigation tracking coordinate system to the medical imaging physical coordinate system, this embodiment of the invention defines rigid body transformation parameters to be solved. The rigid body transformation parameters include a... rotation matrix And one Translation vector .

[0179] The rotation matrix Belongs to a special orthogonal group Satisfying orthogonality constraints And determinant This is used to characterize the attitude rotation deviation between two coordinate systems. The translation vector... Belongs to three-dimensional Euclidean space It is used to characterize the spatial position offset between the origins of two coordinate systems.

[0180] According to the principle of rigid body transformation, the source point set any point in After transformation Represented as:

[0181]

[0182] When solving for the rigid body transformation parameters from the first 3D volume data to the second 3D volume data, the Cholesky method is used to directly calculate the translation parameters when the rigid body transformation only requires translation.

[0183] Step S54: Construct a distance metric model from a point to a tangent plane.

[0184] Unlike Euclidean distance metrics that require complete overlap between points, the metric model constructed in this embodiment of the invention only constrains the transformed source point. Landing at the target point And perpendicular to the normal vector On the tangent plane.

[0185] This invention calculates the transformed source point. To the corresponding number Target tangent plane The normal projection distance. The calculation formula is as follows:

[0186]

[0187] Among them, symbols This represents the dot product operation of vectors. The normal projection distance... Geometrically, this represents the projection component of the position error vector onto the direction of the target normal vector. This is achieved by introducing the tangent plane normal vector. This metric model decouples the normal error from the tangential error, only including the normal error in the registration cost, thus giving the source point the degree of freedom to slide within the target tangential plane.

[0188] Step S55: Establish the weighted least squares global objective function.

[0189] To find the optimal rigid body transformation parameters that minimize the overall registration error of all feature point pairs, this embodiment of the invention constructs a global objective function. The global objective function... Defined as the weighted sum of the squared distances from the points of all valid matching pairs to the tangent plane.

[0190] The specific mathematical expression of the global objective function is as follows:

[0191]

[0192] The symbols are defined as follows: This represents the total number of feature point pairs that effectively participated in the registration. Indicates the first The weight coefficients of the nth feature point pair. In a specific embodiment, the computer processing unit directly calls the weight coefficients of the nth feature point pair obtained in step S41 above. The peak response value of the heatmap for each feature point is used as the weighting coefficient for that point. This setting ensures that anatomical landmarks with higher confidence during feature extraction contribute more to the rigid body transformation parameters during registration optimization, while lower confidence reduces the influence of noise or blurred points. By constructing the above objective function, the rigid body registration problem is transformed into a nonlinear least squares optimization problem. The goal of solving this problem is to find a set of optimal... , so that the objective function value To achieve the minimum, thus realizing the optimal spatial alignment between the surgical navigation coordinate system and the medical imaging coordinate system.

[0193] For rigid body registration problems based on a global objective function, since the rotation matrix belongs to a special orthogonal group... The orthogonality constraint of the objective function makes it nonconvex, and the global optimal solution cannot be obtained directly through analytical methods. Therefore, this embodiment of the invention employs an iterative linearization numerical optimization strategy based on Lie algebra parameterization to solve this nonlinear least squares problem. The specific implementation of the optimization algorithm includes the following steps:

[0194] Step S56: Linearize the rotation variable using Lie algebra or iteratively update it using the conjugate gradient method.

[0195] This invention employs a Gauss-Newton iterative algorithm. First, the rigid body transformation parameters are initialized, for example, by setting the rotation matrix... Let it be the identity matrix, and the translation vector. Set it as the zero vector, or calculate the initial translation using point cloud centroid alignment. In the... In this iteration, the computer processing unit is based on the left perturbation model, assuming the rotation matrix... There is a tiny Lie algebra increment relative to the current estimate. Assume the... The coordinates of the source point at the next iteration are: Apply small rigid body perturbations Updated coordinates The first-order Taylor expansion in the tangent space is approximately expressed as:

[0196]

[0197] Among them, the tiny rigid body disturbance It is The vector, by Rotational Lie algebra vector and translation increment vector Composition; Symbol " "" represents the cross product operation of vectors. Substituting the above linear approximation into the point-to-tangent plane distance formula defined in step S54 above, for the first... For each feature point, the residual equation after linearization is:

[0198]

[0199] In this derivation, the cyclic invariance of the vector mixture product is utilized, and... Equivalent transformation to This allows us to isolate the unknown variables to be solved. ;

[0200] If the conjugate gradient method is used, the processing steps are as follows:

[0201] Calculate the gradient of the objective function: Based on the left perturbation model, calculate the current objective function with respect to the small perturbation vector of the Lie algebra. gradient vector Using the point-to-tangent plane distance formula defined in step S54, for the k-th feature point, its current residual... The source point after transformation The distance to the target tangent plane. According to the chain rule and the properties of the vector mixture product, the gradient contribution generated by a single feature point consists of two parts: [The text abruptly ends here, likely due to an incomplete sentence or missing information.] gradient components: For translation increment vectors gradient components: The gradient contributions of all matched feature points are summed to obtain the current total gradient vector. .

[0202] Determine the conjugate search direction: Calculate the search direction for the current iteration. If it's the first iteration (iter=0), set the search direction to the negative gradient direction: If it is a subsequent iteration (iter>0), calculate the conjugate coefficient. (For example, using the Fletcher-Reeves or Polak-Ribi'ere formula), and utilizing the previous search direction. Construct the current conjugate direction:

[0203] Perform a one-dimensional linear search: along a defined search direction. Perform a one-dimensional search to find an optimal step size scalar. This minimizes the objective function value in that direction. At this point, the optimal Lie algebraic perturbation is determined as... .

[0204] Parameter update: Optimal perturbation amount including rotation and translation components. The incremental transformation matrix is ​​converted to Lie group space through exponential mapping operations and then multiplied on the left by the current rigid body transformation matrix to complete the transformation of the rotation matrix. With translation vector Update.

[0205] Step S57: Construct a weighted linear equation system.

[0206] Based on all The linearized residual equations for each feature point are constructed using the following computer processing unit: The weighted linear equation system. The coefficient matrix It is The matrix whose first... row vector It consists of geometric Jacobian terms relative to the optimization variables, and is specifically defined as follows:

[0207]

[0208] The observation vector It is The vector whose th element Composed of the current residuals, specifically defined as follows:

[0209]

[0210] in, The weighting coefficients are determined in the preceding steps. This is achieved by solving the normal equation. The computer processing unit obtains the optimal perturbation vector that minimizes the sum of squared linearized residuals. In the specific implementation, Cholesky decomposition or Singular Value Decomposition (SVD) algorithms are used to solve this problem. Linear systems are used to ensure numerical stability.

[0211] Step 58: Manifold update of rigid body transformation parameters. After obtaining the optimal perturbation vector, this embodiment of the invention updates the global rigid body transformation parameters, specifically for the translation vector. Perform an addition update in Euclidean space:

[0212]

[0213] For rotation matrix Because of its definition On the manifold, based on the definition of the left perturbation model, this invention example utilizes the Rodriguez formula to transform the rotated Lie algebra vector. Map back to the rotation matrix form and perform a left-multiplication update:

[0214]

[0215] in, This represents the exponential mapping operation from the Lie algebra space to the Lie group space. After the update, the source point coordinates for the next iteration are calculated using the new transformation parameters. Then proceed to the next cycle.

[0216] Step S59: Iterative convergence criteria and result output.

[0217] The optimization algorithm repeats steps S56 to S58 until a preset convergence condition is met. The convergence condition is configured as a combination of the following logic: the norm of the optimal perturbation vector. The value of the rotation matrix is ​​determined when the algorithm converges. The value of the rotation matrix is ​​less than a preset threshold, the relative decrease in the global objective function value is less than a preset tolerance, or the number of iterations reaches a preset upper limit. The translation vector t constitutes the registration matrix required for the surgical navigation system. This registration matrix accurately describes the spatial mapping relationship from the tracking coordinate system where the surgical probe is located to the physical coordinate system where the medical image is located, thus allowing the system to overlay image data onto the patient's actual anatomical structures in real time.

[0218] Based on the same inventive concept as the aforementioned edge feature point extraction and rigid body registration methods, this embodiment of the invention further provides a medical image surgical navigation registration device 100 capable of executing the above-described method flow. This device, through the combination of software functional modules and hardware computing units, achieves automated processing from original image input to final surgical space registration. The specific structure of the medical image surgical navigation registration device 100 includes:

[0219] The dual-branch feature extraction neural network module 110 is configured to load a pre-trained edge feature point extraction network model. Utilizing the parallel computing capabilities of a graphics processing unit (GPU) or a dedicated neural network accelerator, the dual-branch feature extraction neural network module 110 performs three-dimensional convolution operations on the input medical image data. Specifically, the dual-branch feature extraction neural network module 110 internally includes an encoder unit and a decoder unit. The encoder unit is used to extract high-dimensional abstract features of the image layer by layer, and the decoder unit is used to restore spatial resolution and perform multi-task splitting. The data output of the dual-branch feature extraction neural network module 110 is divided into two independent channels: the first channel is connected to the location heatmap output unit, used to output data containing... The device generates a probability response map for each anatomical feature point category; the second channel is connected to the tangential plane pose output unit, which outputs the tangential plane Euler angle parameter tensor corresponding to each voxel position. Through the bi-branch feature extraction neural network module 110, the device converts unstructured medical image data into structured geometric feature descriptors.

[0220] Subvoxel geometry analysis and reconstruction module 120

[0221] The input of the sub-voxel geometry analysis and reconstruction module 120 is connected to the output of the dual-branch feature extraction neural network module 110, used to parse the tensor data output by the network into geometric parameters with physical meaning. The sub-voxel geometry analysis and reconstruction module 120 specifically includes the following processing unit: a local centroid calculation unit 121, used to determine the region of interest based on the peak index of the heatmap, and to perform a weighted average operation using the response values ​​within the local window as weights, thereby calculating sub-voxel level coordinates that break through the image grid resolution limitation. Normal vector reconstruction unit 122 is used to reconstruct the unit normal vector based on the predicted Euler angle parameters using the inverse transformation formula of spherical coordinates. Spatial transformation unit 123 is used to read the affine transformation matrix from the metadata of medical images. The sub-voxel level coordinates are mapped to a new coordinate system using the affine transformation matrix, and the inverse of the affine transformation matrix is ​​used as the final matrix. A covariant vector transformation is performed on the unit normal vector to finally output a set of target tangent planes defined in the patient's physical coordinate system. The semantic matching and data alignment module 130 is used to establish the logical association between the surgical navigation system and medical image data. Specifically, the semantic matching and data alignment module 130 includes the following processing units: a tracer device interface unit 131, used to receive source point sets sent by an external optical tracker or electromagnetic tracker. The source point set contains the coordinates of discrete feature points collected by the doctor using a surgical probe on the surface of the patient's anatomical structure. The index matching unit 132 is used to establish a one-to-one mapping between the points in the source point set and the planar models in the target tangent plane set based on the anatomical category labels of the feature points, constructing a list of matching pairs with definite semantic correspondences. The index matching unit 132 is also configured to filter the matching pairs based on the confidence level of the heatmap, removing feature points with a confidence level lower than a preset threshold.

[0222] The manifold optimization registration solution module 140 is configured to construct and solve a nonlinear least squares optimization problem to achieve high-precision alignment of the spatial coordinate system. Specifically, the manifold optimization registration solution module 140 includes the following processing units: an error metric construction unit 141, used to calculate the normal projection distance from the transformed source point to the target tangent plane. The error metric construction unit 141 uses a point-to-tangent plane distance model, constraining only the deviation in the normal direction, allowing the source point to slide within the tangent plane; and an iterative linearization solution unit 142, used to parameterize the rigid body rotation matrix into a Lie algebra vector based on Lie group and Lie algebra theory, and construct a Gauss-Newton iterative framework. In each iteration cycle, the iterative linearization solution unit 142 calculates the Jacobian matrix and residual vector, solves the linearized normal equation to obtain the optimal perturbation, and updates the global rotation matrix and translation vector through exponential mapping. The output of the manifold optimization registration solution module 140 is used to output the finally converged rigid body transformation matrix. The rigid body transformation matrix is ​​transmitted to the surgical navigation visualization interface to drive the real-time overlay display of medical images and the patient's actual anatomical structure.

[0223] From a hardware perspective, the functional logic of the aforementioned dual-branch feature extraction neural network module 110, sub-voxel geometry analysis and reconstruction module 120, semantic matching and data alignment module 130, and manifold optimization registration solution module 140 is stored in memory 150 and executed by processor 160. Memory 150 stores computer-executable instructions and the weight parameter file of the edge feature point extraction network. Memory 150 includes, but is not limited to, high-speed random access memory (RAM) and non-volatile memory. Processor 160 is connected to memory 150 via a system bus and is configured to execute the computer-executable instructions to implement the functions defined by the aforementioned modules. Processor 160 can be a general-purpose central processing unit, digital signal processor, application-specific integrated circuit, field-programmable gate array, or other programmable logic device. In a specific electronic device embodiment, processor 160 is also connected to input / output interfaces for connecting to medical imaging scanning equipment to acquire image data, connecting to surgical navigation tracking equipment to acquire probe coordinates, and connecting to a display to present the registered navigation screen. Those skilled in the art will understand that the electronic device can be a standalone surgical navigation workstation or an embedded computing unit integrated into a medical robot control cabinet.

[0224] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for edge feature point extraction and registration for rigid body registration of three-dimensional volume data, characterized in that, Includes the following steps: Step S1: Obtain the first three-dimensional volume data and the second three-dimensional volume data to be registered. The first three-dimensional volume data is defined as the source data, and the second three-dimensional volume data is defined as the target data. Step S2: Input the first 3D volume data and the second 3D volume data into the pre-trained UNet network respectively, and output the corresponding edge feature point heatmap and tangent plane pose parameter map; Step S3: For each edge feature point category, analyze the geometric features and construct a tangent plane. For each edge feature point category, use the edge feature point heatmap to determine the edge feature point position, and use the tangent plane pose parameter map to construct the tangent plane normal vector. Step S4: Establish a one-to-one matching relationship between the edge feature points of the first three-dimensional volume data and the second three-dimensional volume data based on the edge feature point categories. The edge feature points in the first three-dimensional volume data are directly matched with the edge feature points of the same category in the second three-dimensional volume data to form matching point pairs. Step S5: Construct an objective function based on the distance from the point to the tangent plane, and solve for the rigid body transformation parameters from the first three-dimensional volume data to the second three-dimensional volume data.

2. The edge feature point extraction and registration method for rigid body registration of three-dimensional volume data according to claim 1, characterized in that, The UNet network is constructed as a fully convolutional neural network structure. The specific process of outputting the corresponding edge feature point heatmap and tangent plane pose parameter map in step S2 includes the following steps: An edge feature point heatmap is generated through a first output branch, which is configured with a three-dimensional convolutional layer and a Sigmoid non-linear activation function layer. Each channel of the edge feature point heatmap corresponds to a class of edge feature points with fixed anatomical semantics. The tangential plane attitude parameter map is generated through the second output branch, which is configured with a three-dimensional convolutional layer, a hyperbolic tangent activation function layer, and a scaling operation. The values ​​of the three channels of the tangential plane attitude parameter map represent the pitch angle, yaw angle, and roll angle in the Euler system, respectively.

3. The edge feature point extraction and registration method for rigid body registration of three-dimensional volume data according to claim 1, characterized in that, The specific process of determining the location of edge feature points using the edge feature point heatmap in step S3 includes the following steps: A full-image maximum search is performed in the corresponding channel of the edge feature point heatmap, and the voxel coordinates with the largest response value are selected as integer-level peak coordinates; A local cube window is defined centered on the integer-level peak coordinates. The values ​​of the predicted heatmap within the local cube window are used as weights to perform a weighted average operation, and the sub-voxel level coordinates are calculated.

4. The edge feature point extraction and registration method for rigid body registration of three-dimensional volume data according to claim 3, characterized in that, Step S3 further includes a process of performing physical space coordinate mapping and restoration on the sub-voxel level coordinates, including the following steps: Read the affine transformation matrix from the metadata of the first three-dimensional volume data or the second three-dimensional volume data. The affine transformation matrix describes the transformation relationship from the image array coordinate system to the DICOM patient coordinate system. The subvoxel-level coordinates are mapped using the affine transformation matrix to obtain the coordinates of the physical feature points defined in the DICOM patient coordinate system.

5. The edge feature point extraction and registration method for rigid body registration of three-dimensional volume data according to claim 1, characterized in that, The specific process of constructing the tangent plane normal vector using the tangent plane attitude parameter map in step S3 includes the following steps: At the coordinates corresponding to the edge feature point positions in the tangential plane attitude parameter diagram, read the corresponding pitch angle, yaw angle, and roll angle. A rotation matrix is ​​constructed according to the preset rotation order of ZYX. The rotation matrix is ​​obtained by multiplying the matrix of the roll angle around the Z-axis, the yaw angle around the Y-axis, and the pitch angle around the X-axis. A preset reference axis vector is selected, and the reference axis vector is transformed using the rotation matrix to calculate the normal vector of the tangent plane.

6. The edge feature point extraction and registration method for rigid body registration of three-dimensional volume data according to claim 5, characterized in that, The process of constructing the tangent plane normal vector using the tangent plane pose parameter map also includes the step of transforming the tangent plane normal vector to the DICOM patient coordinate system: Read the affine transformation matrix of the first three-dimensional volume data or the second three-dimensional volume data, and extract the rotation and scaling matrix from the affine transformation matrix. Both are defined in spatial position in the DICOM patient coordinate system. Calculate the inverse of the rotation-scaling matrix; The physical tangent plane normal vector is obtained by transforming the tangent plane normal vector using the inverse matrix and performing a normalization operation.

7. The edge feature point extraction and registration method for rigid body registration of three-dimensional volume data according to claim 1, characterized in that, The edge feature point category has clear medical semantic attributes and is defined as strictly unique in each three-dimensional volume data; The physical location corresponding to the edge feature point category is located on the surface boundary of human organs or tissues; the process of establishing the matching point pair in step S4 does not perform similarity calculation based on image grayscale or local descriptors.

8. The edge feature point extraction and registration method for rigid body registration of three-dimensional volume data according to claim 1, characterized in that, The specific process of constructing the objective function based on the distance from the point to the tangent plane in step S5 includes: Calculate the distance from the edge feature points of the first 3D volume data after transformation by the rigid body transformation parameters to the normal direction of the tangent plane of the corresponding edge feature points of the second 3D volume data, and sum the squares of the distances to construct an objective function based on the distance from the point to the tangent plane.

9. The edge feature point extraction and registration method for rigid body registration of three-dimensional volume data according to claim 1, characterized in that, The specific process of solving the rigid body transformation parameters from the first three-dimensional volume data to the second three-dimensional volume data in step S5 includes: When a rigid body transformation requires both translation and rotation, a Gauss-Newton iterative algorithm or conjugate gradient method based on Lie algebra parameterization is used to parameterize the rotation matrix into a Lie algebra vector or a residual vector. In each iteration step, a weighted linear equation system is constructed, the Jacobian matrix and residual vector are calculated, and the linearized normal equation is solved to obtain the optimal perturbation amount. The rotation matrix and the translation vector are updated using the exponential mapping operation based on the optimal perturbation amount until the preset convergence condition is met.

10. The edge feature point extraction and registration method for rigid body registration of three-dimensional volume data according to claim 9, characterized in that, The specific process of determining the preset convergence condition in step S5 includes the following steps: Calculate the overall numerical amplitude of the optimal perturbation amount and compare the overall numerical amplitude with a preset minimum error threshold; Simultaneously monitor the number of iterations in the iterative calculation and compare the number of iterations with a preset maximum iteration limit; When the overall numerical amplitude is less than the minimum error threshold, or when the number of iterations reaches the maximum iteration limit, the preset convergence condition is satisfied, and the iterative calculation is terminated.

11. The method for edge feature point extraction and registration for rigid body registration of three-dimensional volume data according to claim 9, characterized in that, The specific process of solving the rigid body transformation parameters from the first three-dimensional volume data to the second three-dimensional volume data in step S5 includes: When a rigid body transformation only requires translation, the Cholesky method can be used to directly calculate the translation parameters.

Images