Medical image registration method based on self-supervised representation learning and hierarchical feature optimization

Abstract

The invention discloses a medical image registration method based on self-supervised representation learning and hierarchical feature optimization, and the method comprises the steps: carrying out the pre-training of a 3D-PVT encoder through a mask auto-encoder strategy, and learning the feature representation of an anatomical structure from a three-dimensional reference image without a label; the method comprises the steps of receiving a reference image and a to-be-registered image of a patient, extracting a multi-level feature pyramid through a pre-trained 3D-PVT encoder, executing coarse-to-fine iterative optimization on the feature pyramid layer by layer on the basis of a mutual information maximization target, selecting an elastic registration mode or a rigid registration mode according to registration requirements in the optimization process, and outputting a registration result. According to the invention, high-precision elastic or rigid registration without manual labeling can be realized.

Description

Technical Field

[0001] This invention belongs to the field of medical image registration technology, specifically relating to a medical image registration method based on self-supervised representation learning and hierarchical feature optimization. Background Technology

[0002] Medical image registration is a key technology in medical image analysis. Its purpose is to spatially align images from different times, devices, or modalities to achieve precise lesion comparison, effective information fusion, and accurate guidance for treatment plans. In the field of radiotherapy, accurate image registration is particularly important. For example, registering planned CT (Computed Tomography) with daily cone-beam CT (CBCT) is fundamental to ensuring accurate radiation dose delivery to the target area while protecting surrounding healthy tissue.

[0003] In clinical practice, registration requirements fall into two categories: rigid registration, primarily used to correct global translation and rotation caused by positioning errors, and commonly used for aligning rigid structures such as bones; and non-rigid (or elastic) registration, used to align soft tissues that deform due to physiological activities such as respiration, organ peristalsis, or changes in filling status. CBCT images, in particular, typically exhibit a lower signal-to-noise ratio, more severe scattering artifacts, and different grayscale distributions compared to planning CT, posing a significant challenge to accurate CT-CBCT registration.

[0004] Existing registration methods mainly include:

[0005] The first category is registration methods based on end-to-end deep learning (represented by VoxelMorph):

[0006] This type of method defines the registration problem as a learning task, using a convolutional neural network (CNN) to construct a function that directly maps input image pairs to a deformation field that aligns the images. The training process is unsupervised, training the network parameters by optimizing a loss function consisting of image similarity loss (such as mean squared error (MSE) or cross-correlation (CC)) and a smoothness regularization term. A significant advantage of this type of method is its extremely fast inference speed, as it transforms the time-consuming optimization process for each image pair in traditional methods into a single forward propagation of the network. Some studies have also explored using auxiliary segmentation labels during training to improve registration accuracy.

[0007] However, this type of method has the following problems:

[0008] 1) Limited feature representation capability: The traditional convolutional neural network (CNN) architecture it adopts has inherent limitations in capturing long-range, global anatomical structural dependencies. The receptive field of the network needs to be gradually expanded by stacking multiple layers of convolutions, which may lead to its insufficient perception of large-scale, non-local complex deformations.

[0009] 2) Direct mapping: The image intensity is directly mapped to the deformation field end-to-end, making the model more sensitive to images outside the training data distribution (such as different artifacts or noise patterns). The optimization process is like a "black box" and has poor interpretability.

[0010] The second category is multimodal registration methods based on anatomical structure perception (represented by CN120070521A):

[0011] This approach aims to address the registration challenge of multimodal images (such as CT and MRI) by enhancing the network's ability to perceive anatomical structures. Its core innovation lies in the design of the "Anatomy-Aware Module (ASPM)," which implicitly embeds traditional image processing operators such as Sobel and Laplacian into convolutional layers, forcing the network to extract first- and second-order differential features of the image (such as edges and gradients). Furthermore, a dedicated edge loss term is added to its loss function to constrain the alignment of anatomical structure edges between images from different modalities.

[0012] However, this type of method has the following problems:

[0013] Reliance on manually designed features: This method relies on predefined, manually designed operators such as Sobel and Laplacian to extract anatomical features. While these fixed operators are interpretable, they may not be optimal choices, as they can only capture simple edge and gradient information and may ignore more complex contextual features such as texture and shape that are crucial for registration. Its performance is limited by these manually designed prior knowledge.

[0014] The third category is multimodal registration methods based on cascaded networks and sequence information (represented by CN114119689B):

[0015] This scheme primarily targets the registration of multimodal MR images. Its key feature is the use of a cascaded network structure, employing two sub-networks (N1 and N2) for coarse-to-fine registration. To address the issue of insufficient training data, this method utilizes the Bloch equation to simulate a large number of distorted training samples. Another crucial feature is the inclusion of sequence information of the MR images (such as echo time and resolution) as an additional input channel, fed into the network along with the images. This aims to help the network cope with differences between different MR sequences.

[0016] However, this type of method has the following problems:

[0017] 1) Fidelity dependent on simulation: The effectiveness of this method is highly dependent on the similarity between simulation data and real clinical data. If the simulation process (such as the Bloch equation) cannot fully reproduce the complex noise and artifact patterns in real clinical images (especially CBCT), then the performance of the model trained on simulation data on real data will be greatly reduced.

[0018] 2) The framework lacks versatility: This method introduces parameters specific to MR sequences as network inputs, which makes its design highly coupled with specific modalities, making it difficult to directly extend or apply to other multimodal registration tasks that do not have such parameters, such as CT-CBCT registration. Summary of the Invention

[0019] To address the aforementioned technical problems, this invention proposes a medical image registration method based on self-supervised representation learning and hierarchical feature optimization.

[0020] To achieve the above objectives, the technical solution of the present invention is as follows:

[0021] In a first aspect, this invention discloses a medical image registration method based on self-supervised representation learning and hierarchical feature optimization, comprising:

[0022] Step S1: Pre-train the 3D-PVT encoder using a masked autoencoder strategy to learn anatomical feature representations from unlabeled 3D reference images;

[0023] Step S2: Receive the patient's baseline image and the image to be registered. Extract a multi-level feature pyramid using a pre-trained 3D-PVT encoder. Based on the mutual information maximization objective, perform iterative optimization layer by layer from coarse to fine on the feature pyramid. During the optimization process, select either an elastic registration mode or a rigid registration mode according to the registration requirements, and output the registration result.

[0024] Based on the above technical solution, the following improvements can be made:

[0025] As a preferred embodiment, step S1 includes:

[0026] Step S1.1: Acquire unlabeled 3D reference images;

[0027] Step S1.2: Preprocess the acquired 3D reference image to unify its voxel spacing;

[0028] Step S1.3: Randomly sample several three-dimensional image blocks from the preprocessed three-dimensional reference image and divide them into non-overlapping voxel cubes, and randomly mask some voxel cubes at a certain ratio.

[0029] Step S1.4: Use a 3D-PVT encoder to process the visible voxel cubes, reconstruct the masked voxel cubes using a decoder, and continuously optimize the 3D-PVT encoder with the goal of minimizing the difference between the reconstructed voxel cubes and the original voxel cubes.

[0030] As a preferred embodiment, step S2 includes:

[0031] Step S2.1: Preprocess the reference image and the image to be registered to unify their voxel spacing so that it meets the input requirements of the 3D-PVT encoder;

[0032] Step S2.2: Using a sliding window strategy, pairs of image patches are extracted from corresponding positions of the preprocessed reference image and the image to be registered in an overlapping manner;

[0033] Step S2.3: Input the paired image blocks into the pre-trained 3D-PVT encoder. The 3D-PVT encoder outputs the feature pyramid of the reference image and the feature pyramid of the image to be registered.

[0034] Step S2.4: Perform feature fusion on the multi-channel reference feature maps output from each level of the feature pyramid of the reference image to obtain the single-channel reference feature map of the corresponding level.

[0035] Feature fusion is performed on the multi-channel feature maps of the image to be registered at each level of the feature pyramid to obtain the single-channel feature map of the corresponding level.

[0036] Step S2.5: Starting from the coarsest layer of the feature pyramid, iterate layer by layer to the finest layer, and find the optimal transformation parameters for each layer. The optimal transformation parameters can maximize the mutual information between the single-channel reference feature map and the transformed single-channel feature map to be registered. During the optimization process, select the flexible registration mode or the rigid registration mode according to the registration requirements, and output the registration result.

[0037] As a preferred embodiment, the mathematical expression for mutual information in step S2.5 is as follows:

[0038]

[0039] The optimization objective function for the i-th layer of the feature pyramid is as follows:

[0040]

[0041] in:

[0042] i is an integer representing the current level index of the feature pyramid;

[0043] T i The transformation parameters to be optimized at the i-th level;

[0044] This is a baseline feature map for a single channel;

[0045] This is a single-channel feature map to be registered;

[0046] o is a transformation composite operator;

[0047] MI(A,B) is a function that calculates the mutual information between two single-channel feature maps A and B.

[0048] p(a,b) is the joint probability distribution of the intensity values ​​a and b at the corresponding voxel positions in the feature map;

[0049] p(a) and p(b) are the marginal probability distributions of intensity values ​​a and b in the feature map, respectively;

[0050] lambda is a scalar hyperparameter;

[0051] For application to transformation T i The smoothness regularization function on the surface.

[0052] As a preferred option, in step S2.5, when the registration requirement selects the flexible registration mode, the optimization process includes:

[0053] Step A1: Transform the parameter T to be optimized i Defined as dense deformation field Φ i , and Φ i :R 3 ->R 3 It maps each voxel coordinate p within the domain to a unique displacement vector d;

[0054] Step A2: Based on the dense deformation field Φ i The feature map to be registered is processed by the displacement mapping of the following formula to obtain the deformed feature map to be registered;

[0055] p ′ =p+Φ i (p);

[0056] Where: p is the feature map to be registered. any point in it;

[0057] p′ is the new coordinate after deformation, calculated by vector addition;

[0058] Step A3: Add regular expression terms Defined as deformation field Φ i L2 norm of spatial gradient:

[0059]

[0060] Where: grad(Φ i (p)): Deformation field Φ i The spatial gradient at voxel p has a norm that measures the degree of local stretching or compression at that point.

[0061] The deformation field Φ is directly adjusted using a numerical optimization algorithm. i The components of the displacement vector are used to maximize the optimization objective function defined in step S2.5;

[0062] Step A4: After optimization and global reconstruction at all levels, the final output is a global dense deformation field Φ covering the entire image size. final .

[0063] As a preferred embodiment, step A4 includes:

[0064] Step A4.1: Create a zero tensor with the same size as the original image to be registered: the globally accumulated deformation field Φ global And the global cumulative weight graph W global ;

[0065] Step A4.2: Traverse all N local deformation fields Φ1,Φ2,…,Φ generated by the sliding window. N};

[0066] For any local deformation field Φ k Its spatial domain in the global image is Ω. k Define a three-dimensional Gaussian weighted window function w k (p local )as follows:

[0067]

[0068] Where, p local Ω k Local coordinates within;

[0069] p center The coordinates of the center of the window;

[0070] σ is the standard deviation controlling the rate of weight decay;

[0071] For Ω k For each global coordinate p within the range, perform an accumulation operation to update the global accumulated deformation field Φ. global And the global cumulative weight graph W global :

[0072] Φ global (p)=Φ global (p)+w k (p local )*Φk (p local );

[0073] W global (p)=W global (p)+w k (p local );

[0074] Step A4.3: After traversing all local deformation fields, normalize the accumulated results to calculate the final global deformation field Φ. final ,

[0075]

[0076] As a preferred option, in step S2.5, when the registration requirement selects the rigid registration mode, the optimization process includes:

[0077] Step B1: Transform the parameter T to be optimized i Defined as a parameter vector θ with 6 degrees of freedom i ;

[0078] θ i ={t x ,t y ,t z ,α,β,γ);

[0079] Among them, t x ,t y ,t z α, β, γ represent the translation along the three coordinate axes, and α, β, γ represent the rotation angles around the three coordinate axes.

[0080] Step B2: Based on the parameter vector θ i Construct the affine transformation matrix M(θ) i The following mapping is used to process the feature map to be registered, resulting in the deformed feature map to be registered;

[0081] p′ h =M(θ) i )*p h ;

[0082] p h =[p x ,p y ,p z ,1] T ;

[0083] M(θ i )=[[R(α,β,γ),[t x ,t y ,t z ] T ],[0,0,0,1]];

[0084] Where: p h Feature map to be registered The homogeneous coordinates of any point in the array;

[0085] p′ h The new coordinates after deformation are calculated by matrix multiplication;

[0086] R(α,β,γ) is a 3x3 rotation matrix generated by the rotation angle;

[0087] Step B3: Adjust θ using a numerical optimization algorithm i The six parameters are used to maximize the optimization objective function defined in step S2.5;

[0088] Step B4: Predict the rigid parameter set {θ1,θ2,...,θ} in blocks using a sliding window strategy. M}, and integrate the rigid parameter set to output the global rigid transformation matrix M. final .

[0089] Secondly, this invention discloses a medical image registration device based on self-supervised representation learning and hierarchical feature optimization, comprising:

[0090] The pre-training module is used to pre-train the 3D-PVT encoder using a masked autoencoder strategy to learn anatomical feature representations from unlabeled 3D reference images.

[0091] The registration module receives the patient's baseline image and the image to be registered. It extracts a multi-level feature pyramid through a pre-trained 3D-PVT encoder and performs iterative optimization layer by layer on the feature pyramid based on the mutual information maximization objective. During the optimization process, it selects either an elastic registration mode or a rigid registration mode according to the registration requirements and outputs the registration result.

[0092] Thirdly, the present invention discloses a computing device, comprising:

[0093] One or more processors;

[0094] Memory;

[0095] And one or more programs, wherein the one or more programs are stored in memory and configured to be executed by one or more processors, and the one or more programs include instructions for any of the above-described medical image registration methods based on self-supervised representation learning and hierarchical feature optimization.

[0096] Fourthly, the present invention discloses a storage medium storing one or more computer-readable programs, the one or more programs including instructions adapted to be loaded by a memory and executed any of the above-described medical image registration methods based on self-supervised representation learning and hierarchical feature optimization.

[0097] This invention discloses a medical image registration method based on self-supervised representation learning and hierarchical feature optimization, which can achieve high-precision, flexible or rigid registration without manual annotation, and has the following beneficial effects:

[0098] First, this invention uses a 3D-PVT encoder and pre-trains it using a masked autoencoder (MAE), resulting in superior feature representation and higher registration robustness.

[0099] Second, the self-supervised learning strategy adopted in this invention can be trained directly on real clinical CT data without relying on any manually drawn segmentation labels or "gold standard" deformation fields, completely eliminating manual labels and physical simulation, and thus having greater versatility and practicality.

[0100] Third, this invention adopts a hybrid framework design, which innovatively decouples deep learning from classical optimization theory, taking into account accuracy, efficiency and interpretability.

[0101] Fourth, the framework of this invention has inherent flexibility. By switching the type of transformation to be solved (high-degree-of-freedom dense deformation field and low-degree-of-freedom 6-parameter rigid transformation) in the optimization step, two functions can be achieved within the same technical framework. Attached Figure Description

[0102] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0103] Figure 1 A flowchart of a medical image registration method provided in an embodiment of the present invention. Detailed Implementation

[0104] The preferred embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0105] The technical solutions of 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.

[0106] The expression “includes” is an “open-ended” expression, which means that there is a corresponding component or step, and should not be interpreted as excluding additional components or steps.

[0107] To achieve the objectives of this invention, in some embodiments of the medical image registration method based on self-supervised representation learning and hierarchical feature optimization, the medical image registration method aims to achieve high-precision, flexible or rigid registration without manual annotation, including the following two major stages.

[0108] Phase 1: Pre-training phase of the self-supervised feature encoder;

[0109] Phase 2: Online image registration based on hierarchical feature maps and mutual information.

[0110] like Figure 1 As shown, the medical image registration method specifically includes:

[0111] Step S101: Pre-train the 3D-PVT encoder using a masked autoencoder strategy to learn anatomical feature representations from unlabeled 3D reference images;

[0112] Step S102: Receive the patient's baseline image (e.g., planned CT image) and the image to be registered (e.g., CBCT image). Extract a multi-level feature pyramid using a pre-trained 3D-PVT encoder. Based on the mutual information maximization objective, perform iterative optimization layer by layer from coarse to fine on the feature pyramid. During the optimization process, select either an elastic registration mode or a rigid registration mode according to the registration requirements, and output the registration result.

[0113] Each of the above steps will be explained in detail below.

[0114] Step S101 includes:

[0115] Step S101.1: Acquire unlabeled 3D reference images;

[0116] Step S101.2: Preprocess the acquired 3D reference image to unify its voxel spacing;

[0117] Step S101.3: Randomly sample several three-dimensional image blocks from the preprocessed three-dimensional reference image and divide them into non-overlapping voxel cubes, and randomly mask some voxel cubes at a certain ratio.

[0118] Step S101.4: Use a 3D-PVT encoder to process the visible voxel cubes, reconstruct the masked voxel cubes using a decoder, and continuously optimize the 3D-PVT encoder with the goal of minimizing the difference between the reconstructed voxel cubes and the original voxel cubes.

[0119] In a specific embodiment, the aforementioned three-dimensional reference image may be, but is not limited to, a clinical three-dimensional CT image.

[0120] To eliminate the interference of voxel spacing differences caused by different equipment and scanning parameters on the model's spatial feature learning, all CT images need to be preprocessed. Isotropic resampling is performed using methods such as trilinear interpolation to unify the voxel spacing of all data to 1mm x 1mm x 1mm. This step is crucial to ensuring that the model learns a consistent and accurate spatial representation of anatomical structures.

[0121] This invention employs a 3D Pyramid Vision Transformer (3D-PVT) as the backbone network for feature extraction. The advantage of this network architecture lies in its ability to naturally output a feature map pyramid (e.g., a four-stage feature map) containing multiple spatial resolutions and semantic levels during a single forward propagation, providing direct input for subsequent hierarchical registration strategies.

[0122] Step S101.3 adopts a self-supervised learning strategy of Masked Autoencoder (MAE).

[0123] Specifically, three-dimensional image blocks (e.g., 96x96x96 voxels) are randomly sampled from preprocessed CT data. After the image blocks are segmented into a series of non-overlapping voxel cubes, most of the voxel cubes are randomly masked at a high percentage (e.g., 75%).

[0124] The 3D-PVT encoder only receives and processes the visible, unmasked voxel cubes and outputs their high-dimensional feature representations. A lightweight decoder receives the encoder's output and attempts to reconstruct the original masked voxel cubes based on this limited feature information.

[0125] It is worth noting that the aforementioned pre-training method does not rely on any manually labeled segmentation or registration tags. Instead, it drives the model to learn the inherent spatial relationships and semantic information of anatomical structures by allowing the model to "fill in the blanks" on its own. Compared with existing techniques that rely on labeled data or simply use autoencoders, MAE enables the model to learn more robust and higher-level feature representations, greatly improving the model's generalization ability.

[0126] Specifically, a three-dimensional image block (e.g., 96x96x96 voxels) is randomly sampled from the preprocessed CT data and denoted as x.

[0127] The image patch x is divided into a set P = {p1, p2, ..., p...} of non-overlapping voxel cubes. N}

[0128] Then, this set is randomly sampled at a relatively high percentage (e.g., 90%) to obtain the visible subset P. vis and the masked subset P mask .

[0129] The 3D-PVT encoder E only receives and processes the visible subset P. vis And map it to a compact high-dimensional latent representation z vis =E(P vis ).

[0130] A lightweight decoder D receives the encoder's output z_vis and, combined with the position information of the masked cube, attempts to reconstruct the original masked voxel cube content P′. mask =D(z) vis ).

[0131] The optimization objective of this process is to minimize the difference between the reconstructed voxel cube and the original voxel cube. Mean Squared Error (MSE) is used as the loss function L. MAE Its mathematical expression is:

[0132]

[0133] Among them, |P mask | represents the number of masked voxel cubes, p i It is the original shielded voxel cube, p ′ i It is the corresponding voxel cube reconstructed by the decoder.

[0134] After pre-training, the decoder is discarded, and only the trained 3D-PVT encoder is retained. The weights of this encoder are "frozen" to serve as a fixed, efficient feature extraction module for the next stage of online registration.

[0135] This invention discloses a self-supervised medical image feature learning method based on masked autoencoder (MAE) and 3D pyramid visual transformer (3D-PVT). The method uses MAE to pre-train the ViT architecture to learn high-quality anatomical feature representations with global receptive fields from unlabeled 3D medical images.

[0136] Furthermore, this invention adopts a hybrid registration framework that decouples feature extraction and registration optimization, dividing the process into two independent stages: feature extraction is performed using a pre-trained deep network, and then registration is performed on the feature map using a numerical optimization algorithm based on the classic similarity metric (mutual information), instead of using an end-to-end learning paradigm.

[0137] Step S102 includes:

[0138] Step S102.1: Preprocess the reference image and the image to be registered to unify their voxel spacing so that they meet the input requirements of the 3D-PVT encoder;

[0139] Step S102.2: Using a sliding window strategy, pairs of image patches are extracted from corresponding positions of the preprocessed reference image and the image to be registered in an overlapping manner;

[0140] Step S102.3: Input paired image blocks into the pre-trained 3D-PVT encoder. The 3D-PVT encoder outputs the feature pyramid of the reference image and the feature pyramid of the image to be registered.

[0141] Step S102.4: Perform feature fusion on the multi-channel reference feature maps output from each level of the feature pyramid of the reference image to obtain the single-channel reference feature map of the corresponding level.

[0142] Feature fusion is performed on the multi-channel feature maps of the image to be registered at each level of the feature pyramid to obtain the single-channel feature map of the corresponding level.

[0143] Step S102.5: Starting from the coarsest layer of the feature pyramid, iterate layer by layer to the finest layer, and find the optimal transformation parameters for each layer. The optimal transformation parameters can maximize the mutual information between the single-channel reference feature map and the transformed single-channel feature map to be registered. During the optimization process, select the flexible registration mode or the rigid registration mode according to the registration requirements, and output the registration result.

[0144] This invention discloses a hierarchical optimization registration strategy based on multi-scale deep feature maps. It utilizes the feature pyramid output by PVT to perform iterative optimization from the coarsest feature level upwards layer by layer, and solves the transformation that maximizes mutual information at each layer.

[0145] Specifically, for complete, large-sized input images (e.g., 512x512x128 voxels), a sliding window strategy is used for processing.

[0146] First, isotropic resampling is performed using methods such as trilinear interpolation to match the spatial resolution of the training dataset, ensuring that the model can correctly extract feature information.

[0147] Then, paired image patches (96x96x96 voxels) are extracted from corresponding positions in the reference image and the image to be registered in an overlapping manner (e.g., with a step size of 48 voxels and an overlap rate of 50%).

[0148] This invention employs an iterative framework from coarse to fine, which optimizes the transformation relationship layer by layer on the feature pyramid output by the pre-trained encoder.

[0149] Step S102.4 Dimensionality reduction of the feature channels. The encoder outputs the baseline feature map at each pyramid level i. and the feature map to be registered All are multi-channel feature maps, with dimensions of (height * width * depth * number of channels C). To simplify the joint probability distribution estimation in subsequent mutual information calculation and to extract the most critical fusion features, a weighted average is performed on the input C channel feature maps, thereby effectively fusing the multi-channel feature maps into a single-channel feature map.

[0150] This step yields the baseline feature map for a single channel. and the feature map to be registered Subsequent optimizations were all performed on these single-channel feature maps.

[0151] In the subsequent registration optimization process of mutual information, directly estimating the joint probability distribution on the high-dimensional multi-channel feature map will significantly increase the computational complexity and is easily affected by the "curse of dimensionality," leading to inaccurate statistical estimation. Therefore, this invention uses a channel-weighted average fusion strategy to convert the multi-channel feature map output by the encoder into a single-channel feature map. This operation has the following advantages:

[0152] 1) It can effectively preserve key spatial structural information;

[0153] 2) Reduce computational complexity to make mutual information estimation feasible in real time with 3D data;

[0154] 3) Reduce the interference of random noise on the statistical distribution and improve matching stability.

[0155] 4) The weighting method can be set based on prior knowledge, or a learnable attention mechanism can be adopted in future versions to enhance the adaptive ability.

[0156] This invention discloses a feature channel dimensionality reduction method for medical image registration. Before calculating mutual information, a channel-weighted averaging operation (such as using 1x1x1 convolution) is used to fuse multi-channel depth feature maps into a single-channel feature map, thereby simplifying and stabilizing the subsequent optimization process.

[0157] The optimization process of this invention starts from the coarsest layer of the feature pyramid and iterates to the finest layer. The goal of each layer of the feature pyramid is to solve for an optimal transformation T. i This transformation can maximize the single-channel fixed feature map. Compared with the transformed single-channel motion feature map Mutual information between them.

[0158] Mutual information (MI) is a measure of the statistical dependence between two random variables (in this case, the intensity values ​​of the feature maps), and its mathematical expression is as follows:

[0159]

[0160] To ensure the physical realism of the transformation and prevent unreasonable tearing or folding, a smoothness regularization term L needs to be introduced into the optimization objective. smooth .

[0161] The optimization objective function for the i-th layer of the feature pyramid is as follows:

[0162]

[0163] in:

[0164] i is an integer representing the current level index of the feature pyramid;

[0165] T i The transformation parameters to be optimized at the i-th level can be either elastic deformation fields or rigid transformation parameters, depending on the registration mode.

[0166] This is a baseline feature map for a single channel;

[0167] This is a single-channel feature map to be registered;

[0168] o is the transformation composition operator, indicating that the transformation T will be performed... i Applied to feature maps On the spatial coordinates, to generate a new distorted feature map;

[0169] MI(A,B) is a function that calculates the mutual information between two single-channel feature maps A and B.

[0170] p(a,b) is the joint probability distribution of the intensity values ​​a and b at the corresponding voxel positions in the feature map;

[0171] p(a) and p(b) are the marginal probability distributions of intensity values ​​a and b in the feature map, respectively;

[0172] lambda is a scalar hyperparameter used as a regularization weight to balance registration accuracy (driven by the MI term) and transformation smoothness (driven by L). smooth The relationship between constraints;

[0173] For application to transformation T i The smoothness regularization function on the surface.

[0174] The framework of this invention defines different transformations T to be solved during the optimization process. i It can flexibly support both flexible and rigid registration modes to adapt to different clinical needs.

[0175] Furthermore, in step S102.5, when the registration requirement selects the Elastic Registration mode, the optimization process includes:

[0176] Step A1: Transform the parameter T to be optimized i Defined as dense deformation field Φ i The deformation field is a three-dimensional vector field, and Φ i :R 3 ->R 3 It maps each voxel coordinate p within the domain to a unique displacement vector d;

[0177] Step A2: Based on the dense deformation field Φ i The feature map to be registered is processed by the displacement mapping of the following formula to obtain the deformed feature map to be registered;

[0178] p ′ =p+Φ i (p);

[0179] Where: p is the feature map to be registered. any point in it;

[0180] p′ is the new coordinate after deformation, calculated by vector addition;

[0181] Step A3: Add regular expression terms Defined as deformation field Φ i L2 norm of spatial gradient:

[0182]

[0183] Where: grad(Φ i (p)): Deformation field Φ i The spatial gradient at voxel p has a norm that measures the degree of local stretching or compression at that point.

[0184] Numerical optimization algorithms (such as L-BFGS) can be used to directly adjust the deformation field Φ. i The components of the displacement vector are used to maximize the optimization objective function defined in step S102.5;

[0185] Step A4: After optimization and global reconstruction at all levels, the final output is a global dense deformation field Φ covering the entire image size. final .

[0186] Step A4 above specifically includes:

[0187] Step A4.1: Create a zero tensor with the same size as the original image to be registered:

[0188] Global cumulative deformation field Φ global And the global cumulative weight graph W global ;

[0189] Step A4.2: Traverse all N local deformation fields Φ1,Φ2,…,Φ generated by the sliding window. N};

[0190] For any local deformation field Φ k (Its spatial domain in the global image is Ω) k ), perform the following operations:

[0191] 1) Define a three-dimensional Gaussian weighted window function w with the same size as the window. k (p local This function has a maximum value (e.g., 1) at the center of the window and decays smoothly towards the edges. Its expression is:

[0192]

[0193] Where, p local Ω k Local coordinates within;

[0194] p centerThe coordinates of the center of the window;

[0195] σ is the standard deviation controlling the rate of weight decay;

[0196] 2) For Ω k For each global coordinate p within the range, perform an accumulation operation to update the global accumulated deformation field Φ. global And the global cumulative weight graph W global :

[0197] Φ global (p)=Φ global (p)+w k (p local )*Φ k (p local );

[0198] W global (p)=W global (p)+w k (p local );

[0199] Step A4.3: After traversing all local deformation fields, normalize the accumulated results to calculate the final global deformation field Φ. final For each coordinate p in the global domain Ω:

[0200]

[0201] It is worth noting that the element-by-element division operation here only applies to W. global Regions with p>0 are processed. This weighted average fusion method effectively suppresses boundary artifacts that may be caused by block stitching, ensuring the spatial smoothness and continuity of the final output global deformation field.

[0202] The above steps aim to reduce the large number of local, overlapping deformation fields Φ generated in the sliding window strategy. k Seamlessly merged into a smooth global deformation field Φ covering the entire image domain Ω. final .

[0203] This invention discloses a sliding window inference and weighted fusion method for global images. The method applies a model trained on image patches to a full-size image through overlapping sliding windows, and seamlessly fuses the generated local transformations through Gaussian weighted averaging and other methods to obtain a globally consistent final registration result.

[0204] Furthermore, in step S102.5, when the registration requirement selects the rigid registration mode, the optimization process includes:

[0205] Step B1: Transform the parameter T to be optimizedi Defined as a low-dimensional parameter vector θ i This vector contains rigid transformation parameters with 6 degrees of freedom;

[0206] θ i =(t x ,t y ,t z ,α,β,γ);

[0207] Among them, t x ,t y ,t z α, β, γ represent the translation along the three coordinate axes, and α, β, γ represent the rotation angles (e.g., Euler angles) around the three coordinate axes.

[0208] Step B2: Based on the parameter vector θ i Construct a 4x4 affine transformation matrix M(θ) i The following mapping is used to process the feature map to be registered, resulting in the deformed feature map to be registered;

[0209] p′ h =M(θ) i )*p h ;

[0210] p h =[p x ,p y ,p z ,1] T ;

[0211] M(θ i )=[[R(α,β,γ),[t x ,t y ,t z ] T ],[0,0,0,1]];

[0212] in:

[0213] p h Feature map to be registered The homogeneous coordinates of any point in the array;

[0214] p′ h The new coordinates after deformation are calculated by matrix multiplication;

[0215] R(α,β,γ) is a 3x3 rotation matrix generated by the rotation angle;

[0216] Step B3: Adjust θ using a numerical optimization algorithm i The six parameters are used to maximize the optimization objective function defined in step S102.5;

[0217] Step B4: Predict the rigid parameter set {θ1,θ2,...,θ} in blocks using a sliding window strategy. M}, and integrate the rigid parameter set to output the global rigid transformation matrix M. final .

[0218] In sliding window inference mode, each image patch predicts an independent set of stiffness parameters. The final global stiffness transformation matrix M final It can be determined by robust averaging of the parameters predicted for all blocks or by a consensus-based voting algorithm (such as RANSAC).

[0219] Compared to direct global rigidity estimation, this method has the following advantages:

[0220] 1) Stronger local anti-interference capability, which can avoid excessive interference to the overall registration result due to a single area;

[0221] 2) The RANSAC strategy can effectively identify inconsistent or outlier blocks, eliminate abnormal estimates, and improve global stability.

[0222] 3) This method is more adaptable if there are differences in movement between organs or degradation in some areas in the registration scenario.

[0223] This invention discloses a unified technical solution that supports dual-mode registration. Under the same hierarchical optimization framework, elastic registration or rigid registration can be flexibly achieved by selecting different optimization objectives (dense deformation field or 6-DOF rigid transformation parameters).

[0224] In other embodiments, the present invention discloses a medical image registration device based on self-supervised representation learning and hierarchical feature optimization, comprising:

[0225] The pre-training module is used to pre-train the 3D-PVT encoder using a masked autoencoder strategy to learn anatomical feature representations from unlabeled 3D reference images.

[0226] The registration module receives the patient's baseline image and the image to be registered. It extracts a multi-level feature pyramid through a pre-trained 3D-PVT encoder and performs iterative optimization layer by layer on the feature pyramid based on the mutual information maximization objective. During the optimization process, it selects either an elastic registration mode or a rigid registration mode according to the registration requirements and outputs the registration result.

[0227] Furthermore, the pre-training module includes:

[0228] Acquisition unit, used to acquire label-free 3D reference images;

[0229] The pre-training preprocessing unit is used to preprocess the acquired 3D reference images to unify their voxel spacing;

[0230] The segmentation unit is used to randomly sample several three-dimensional image blocks from the preprocessed three-dimensional reference image and segment them into non-overlapping voxel cubes, and randomly mask some voxel cubes at a certain ratio.

[0231] The pre-trained optimization unit is used to process visible voxel cubes using the 3D-PVT encoder, reconstruct masked voxel cubes through the decoder, and continuously optimize the 3D-PVT encoder with the goal of minimizing the difference between the reconstructed voxel cube and the original voxel cube.

[0232] Furthermore, the registration module includes:

[0233] The registration preprocessing unit is used to preprocess the reference image and the image to be registered to unify their voxel spacing so that they meet the input requirements of the 3D-PVT encoder.

[0234] The sliding window extraction unit is used to extract paired image patches from corresponding positions of the preprocessed reference image and the image to be registered in an overlapping manner using a sliding window strategy.

[0235] The encoding unit is used to input pairs of image blocks into a pre-trained 3D-PVT encoder, and the 3D-PVT encoder outputs the feature pyramid of the reference image and the feature pyramid of the image to be registered.

[0236] The dimension reduction unit is used to perform feature fusion on the multi-channel baseline feature maps output by each level of the feature pyramid of the baseline image to obtain the single-channel baseline feature map of the corresponding level.

[0237] Feature fusion is performed on the multi-channel feature maps of the image to be registered at each level of the feature pyramid to obtain the single-channel feature map of the corresponding level.

[0238] The registration optimization unit is used to iterate from the coarsest layer of the feature pyramid to the finest layer, and find the optimal transformation parameters for each layer. The optimal transformation parameters can maximize the mutual information between the single-channel reference feature map and the transformed single-channel feature map to be registered. During the optimization process, the unit selects the flexible registration mode or the rigid registration mode according to the registration requirements and outputs the registration result.

[0239] Furthermore, the mathematical expression for mutual information in the registration optimization unit is as follows:

[0240]

[0241] The optimization objective function for the i-th layer of the feature pyramid is as follows:

[0242]

[0243] in:

[0244] i is an integer representing the current level index of the feature pyramid;

[0245] T i The transformation parameters to be optimized at the i-th level;

[0246] This is a baseline feature map for a single channel;

[0247] This is a single-channel feature map to be registered;

[0248] o is a transformation composite operator;

[0249] MI(A,B) is a function that calculates the mutual information between two single-channel feature maps A and B.

[0250] p(a,b) is the joint probability distribution of the intensity values ​​a and b at the corresponding voxel positions in the feature map;

[0251] p(a) and p(b) are the marginal probability distributions of intensity values ​​a and b in the feature map, respectively;

[0252] lambda is a scalar hyperparameter;

[0253] For application to transformation T i The smoothness regularization function on the surface.

[0254] Furthermore, the registration optimization unit includes: an elastic registration optimization unit, which is used to perform optimization when the registration requirement selects the elastic registration mode;

[0255] Specifically, it includes:

[0256] The elastic parameter definition unit is used to define the transformation parameter T to be optimized. i Defined as dense deformation field Φ i , and Φ i :R 3 ->R 3 It maps each voxel coordinate p within the domain to a unique displacement vector d;

[0257] Elastic processing unit for use based on dense deformation field Φ i The feature map to be registered is processed by the displacement mapping of the following formula to obtain the deformed feature map to be registered;

[0258] p′=p+Φ i (p);

[0259] Where: p is the feature map to be registered. any point in it;

[0260] p′ is the new coordinate after deformation, calculated by vector addition;

[0261] Elastic optimization unit, used to adjust regularization terms Defined as deformation field Φ i L2 norm of spatial gradient:

[0262]

[0263] Where: grad(Φ i (p)): Deformation field Φ i The spatial gradient at voxel p has a norm that measures the degree of local stretching or compression at that point.

[0264] The deformation field Φ is directly adjusted using a numerical optimization algorithm. i The components of the displacement vector are used to maximize the optimization objective function defined in the registration optimization unit;

[0265] The elastic output unit is used to output a globally dense deformation field Φ covering the entire image size after all levels of optimization and global reconstruction. final .

[0266] Furthermore, the elastic output unit includes:

[0267] Initialization unit, used to create a zero tensor with the same size as the original image to be registered: global accumulated deformation field Φ global And the global cumulative weight graph W global ;

[0268] Accumulation unit, used to traverse all N local deformation fields {Φ1,Φ2,…,Φ} generated by the sliding window. N};

[0269] For any local deformation field Φ k Its spatial domain in the global image is Ω. k Define a three-dimensional Gaussian weighted window function w k (p local )as follows:

[0270]

[0271] Where, p local Ω k Local coordinates within;

[0272] p center The coordinates of the center of the window;

[0273] σ is the standard deviation controlling the rate of weight decay;

[0274] For Ω k For each global coordinate p within the range, perform an accumulation operation to update the global accumulated deformation field Φ. global And the global cumulative weight graph W global :

[0275] Φ global (p)=Φ global (p)+w k (p local )*Φ k (p local );

[0276] W global (p)=W global (p)+w k (p local );

[0277] The normalization unit is used to normalize the accumulated results after traversing all local deformation fields in order to calculate the final global deformation field Φ. final ,

[0278]

[0279] Furthermore, the registration optimization unit includes a rigid registration optimization unit, which is used to perform optimization when the registration requirement selects the rigid registration mode;

[0280] Specifically, it includes:

[0281] A rigid parameter definition unit is used to define the transformation parameter T to be optimized. i Defined as a parameter vector θ with 6 degrees of freedom i ;

[0282] θ i =(t x ,t y ,t z ,α,β,γ);

[0283] Among them, t x ,t y ,t z α, β, γ represent the translation along the three coordinate axes, and α, β, γ represent the rotation angles around the three coordinate axes.

[0284] Rigid processing unit, used to process according to parameter vector θ i Construct the affine transformation matrix M(θ) i The following mapping is used to process the feature map to be registered, resulting in the deformed feature map to be registered;

[0285] p′h =M(θ) i )*p h ;

[0286] p h =[p x ,p y ,p z ,1] T ;

[0287] M(θ i )=[[R(α,β,γ),[t x ,t y ,t z ] T ],[0,0,0,1]];

[0288] in:

[0289] p h Feature map to be registered The homogeneous coordinates of any point in the array;

[0290] p′ h The new coordinates after deformation are calculated by matrix multiplication;

[0291] R(α,β,γ) is a 3x3 rotation matrix generated by the rotation angle;

[0292] Rigid optimization elements are used to adjust θ using numerical optimization algorithms. i The six parameters are used to maximize the objective function defined in the registration optimization unit;

[0293] A rigid output unit is used to predict the rigid parameter set {θ1,θ2,...,θ} in blocks using a sliding window strategy. M}, and integrate the rigid parameter set to output the global rigid transformation matrix M. final .

[0294] Furthermore, it should be noted that the medical image registration device based on self-supervised representation learning and hierarchical feature optimization provided in the above embodiments is only illustrated by the division of the above functional modules when performing image registration. In actual applications, the above functions can be assigned to different functional modules as needed, that is, the internal structure of the medical image registration device based on self-supervised representation learning and hierarchical feature optimization can be divided into different functional modules to complete all or part of the functions described above.

[0295] Furthermore, the medical image registration device based on self-supervised representation learning and hierarchical feature optimization provided in the above embodiments and the medical image registration method based on self-supervised representation learning and hierarchical feature optimization belong to the same concept. The specific implementation process can be found in the method embodiments, and will not be repeated here.

[0296] In other embodiments, the present invention discloses a computing device comprising:

[0297] One or more processors;

[0298] Memory;

[0299] And one or more programs, wherein the one or more programs are stored in memory and configured to be executed by one or more processors, and the one or more programs include instructions for any of the above-described medical image registration methods based on self-supervised representation learning and hierarchical feature optimization.

[0300] In other embodiments, the present invention discloses a storage medium storing one or more computer-readable programs, the one or more programs including instructions adapted to be loaded by memory and executed any of the above-described medical image registration methods based on self-supervised representation learning and hierarchical feature optimization.

[0301] This invention discloses a medical image registration method based on self-supervised representation learning and hierarchical feature optimization, which has the following beneficial effects:

[0302] First, it has superior feature representation and higher registration robustness.

[0303] This invention employs a backbone network based on Visual Transformer (3D-PVT) and pre-trains it using a Masked Autoencoder (MAE). Compared to the traditional CNN architecture used in existing technologies such as VoxelMorph, Transformer can more effectively capture the global dependencies and long-range contextual relationships of images. The self-supervised pre-training method of MAE enables the model to learn deeper and more complex anatomical representations from unlabeled data than the manually designed Sobel and Laplacian operators in existing technology CN120070521A.

[0304] Therefore, the features extracted by this invention are more robust to noise, artifacts and large-scale deformation, providing higher quality and more reliable input for subsequent registration optimization, thereby improving the final registration accuracy.

[0305] Second, it completely eliminates manual labeling and physical simulation, making it more versatile and practical.

[0306] The self-supervised learning strategy employed in this invention enables direct training on real clinical CT data, without relying on any manually drawn segmentation labels or "gold standard" deformation fields. This fundamentally differs from supervised learning methods that require large amounts of labeled data. Furthermore, it avoids the data simulation method based on Bloch equations found in existing technology CN114119689B.

[0307] This invention significantly reduces the cost and time required for data preparation and avoids performance degradation caused by inconsistencies in the distribution of simulation data and real clinical data (especially CBCT with complex artifacts). Furthermore, its independence from specific modality sequence information makes its framework more versatile and directly applicable to various registration tasks such as CT-CBCT.

[0308] Third, the hybrid framework design balances accuracy, efficiency, and interpretability.

[0309] This invention innovatively decouples deep learning from classical optimization theory. Unlike end-to-end methods such as VoxelMorph, which directly use the network to predict the deformation field, it utilizes a deep network as a powerful "feature perceptron" and then feeds the high-quality output feature map to a traditional optimization algorithm based on mutual information for precise alignment.

[0310] This hybrid design, on the one hand, utilizes deep features to avoid the problem of traditional methods getting trapped in local optima when optimizing on the original pixels; on the other hand, it retains the advantages of mutual information optimization, which has a clear objective and a solid mathematical foundation, making it more interpretable than "black box" end-to-end methods. While maintaining high accuracy, it improves computational efficiency by optimizing in the feature space rather than the pixel space.

[0311] Fourth, it has a high degree of functional integration and a highly flexible framework.

[0312] The framework of this invention has inherent flexibility. By switching the type of transformation to be solved (high-degree-of-freedom dense deformation field vs. low-degree-of-freedom 6-parameter rigid transformation) in the optimization step, two functions can be achieved within the same technical framework.

[0313] Users can choose to perform high-precision flexible registration to align soft tissues or fast, robust rigid registration to correct positioning, depending on clinical needs. This dual-purpose capability makes it a more comprehensive and practical clinical solution than single-function registration tools.

[0314] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the present invention. Various changes and modifications can be made to the present invention without departing from its spirit and scope. All such changes and modifications fall within the scope of the present invention as claimed, which is defined by the appended claims and their equivalents.

Claims

1. A medical image registration method based on self-supervised representation learning and hierarchical feature optimization, characterized in that, The method comprises the following steps: Step S1: pre-training a 3D-PVT encoder using a mask autoencoder strategy to learn anatomical structure feature representation from unlabeled three-dimensional reference images; Step S2: receiving a reference image and a to-be-registered image of a patient, extracting a multi-level feature pyramid through the pre-trained 3D-PVT encoder, and performing coarse-to-fine iterative optimization on each level of the feature pyramid based on a mutual information maximization target, and selecting an elastic registration mode or a rigid registration mode according to the registration requirement in the optimization process, and outputting a registration result.

2. The medical image registration method of claim 1, wherein, The step S1 comprises: Step S1.1: collecting unlabeled three-dimensional reference images; Step S1.2: preprocessing the collected three-dimensional reference images to unify the voxel spacing thereof; Step S1.3: randomly sampling a plurality of three-dimensional image blocks from the preprocessed three-dimensional reference images, and dividing the three-dimensional image blocks into non-overlapping voxel cubes, and randomly shielding a portion of the voxel cubes at a certain proportion; Step S1.4: processing the visible voxel cubes using the 3D-PVT encoder, reconstructing the shielded voxel cubes through the decoder, and continuously optimizing the 3D-PVT encoder by minimizing the difference between the reconstructed voxel cubes and the original voxel cubes.

3. The medical image registration method of claim 1, wherein, The step S2 comprises: Step S2.1: preprocessing the reference image and the to-be-registered image to unify the voxel spacing thereof, so as to meet the input requirements of the 3D-PVT encoder; Step S2.2: adopting a sliding window strategy to extract paired image blocks from the corresponding positions of the preprocessed reference image and the to-be-registered image in an overlapping manner; Step S2.3: inputting the paired image blocks into the pre-trained 3D-PVT encoder, and outputting a feature pyramid of the reference image and a feature pyramid of the to-be-registered image by the 3D-PVT encoder; Step S2.4: performing feature fusion on the multi-channel reference feature maps output by each level of the feature pyramid of the reference image to obtain single-channel reference feature maps of the corresponding levels; performing feature fusion on the multi-channel to-be-registered feature maps output by each level of the feature pyramid of the to-be-registered image to obtain single-channel to-be-registered feature maps of the corresponding levels; Step S2.5: starting from the coarsest layer of the feature pyramid and iteratively to the finest layer, finding the optimal transformation parameter of each layer, which can maximize the mutual information between the single-channel reference feature maps and the single-channel to-be-registered feature maps after transformation, and selecting an elastic registration mode or a rigid registration mode according to the registration requirement in the optimization process, and outputting a registration result.

4. The medical image registration method according to claim 3, characterized in that, In the step S2.5, the mathematical expression of the mutual information is as follows: The optimization objective function of the i-th layer of the feature pyramid is as follows: Wherein: i is an integer, representing the current level index of the feature pyramid; T i to be optimized on the i-th level; reference feature map is a single-channel reference feature map; is a single-channel feature map to be registered; o is a transformation composite operator; MI(A,B) is a function for calculating the mutual information between two single-channel feature maps A and B; p(a,b) is the joint probability distribution of the intensity values a and b of the corresponding voxel positions in the feature map; p(a), p(b) are the marginal probability distributions of the intensity values a and b in the feature map, respectively; lambda is a scalar hyperparameter; For the smoothness regularizer function applied to the transform T i on the transform T 5. The medical image registration method according to claim 4, characterized in that, In the step S2.5, when the registration requirement selects the elastic registration mode, the optimization process comprises: Step A1: The transformation parameters T to be optimized are defined i defined as a dense deformation field Φ i , and Φ i : R 3 -> R 3 maps each voxel coordinate p in the domain to a unique displacement vector d; Step A2: based on dense deformation field Φ i The to-be-registered feature map is processed by displacement mapping of the following formula to obtain a deformed to-be-registered feature map. p' = p + Φ i (p); wherein: p is any point in the feature map to be registered ​ p' is the new coordinate after deformation calculated by vector addition; Step A3: The regularizer term is defined as the deformation field Φ i L2 norm of spatial gradient: where: grad(Φ i (p)) : deformation field Φ i Spatial gradient at voxel p, whose norm measures the local stretching or compression at that point; using a numerical optimization algorithm to directly adjust the deformation field Φ i components of the displacement vector to maximize the optimization objective function defined in step S2.5; Step A4: After optimization and global reconstruction of all levels, a global dense deformation field Φ covering the full image size is finally output final .

6. The medical image registration method of claim 5, wherein, The step A4 comprises: Step A4.1: Create a zero tensor of the same size as the original image to be registered: global accumulated deformation field Φ global and global accumulated weight map W global ; Step A4.2: Iterate over all N local deformation fields Φ1, Φ2,..., ΦNresulting from the sliding window N}; For any one local deformation field Φ k , which is in spatial domain Ω k in the global image, define a three-dimensional Gaussian weight window function w k (p local ) as follows: where p local is a local coordinate within Ω k ; p center is the center coordinate of the window; σ is the standard deviation of the control weight decay rate; For each global coordinate p in Ω k perform an accumulation operation, updating the global accumulated deformation field Φ global and the global accumulated weight map W global : Φ global (p) = Φ global (p) + w k (p local ) * Φ k (p local ); W global (p) = W global (p) + w k (p local ); Step A4.3: After all local deformation fields have been traversed, the accumulated results are normalized to compute the final global deformation field Φ final , W global (p)>0.

7. The medical image registration method of claim 4, wherein, In the step S2.5, when the registration requirement selects the rigid registration mode, the optimization process comprises: Step B1 : The transformation parameters T to be optimized are defined as a parameter vector θ i defined as a 6-DOF parameter vector θ i ; θ i = (t x , t y , t z , a, b, g); where t x ,t y ,t z represent the translation along the three coordinate axes, and a, b, g represent the rotation angles around the three coordinate axes. Step B2: according to the parameter vector θ i , construct the affine transformation matrix M(θ i ), process the to-be-registered feature map by a mapping of the following formula to obtain a deformed to-be-registered feature map; p' h = M(θ i )*p h ; p h = [p x , p y , p z , 1] T ; M(θ i ) = [[R(α,β,γ),[t x ,t y ,t z ] T ],[0,0,0,1]]; wherein: p h is the homogeneous coordinate of any point in the feature map to be registered; is the homogeneous coordinate of any point in the feature map to be registered; p′ h new coordinates after deformation calculated by matrix multiplication; R(α,β,γ) is a 3x3 rotation matrix generated by the rotation angle; Step B3: Adjust the 6 parameters in Θ using a numerical optimization algorithm to maximize the optimization objective function defined in step S2.

5. i Step B3: Adjust the 6 parameters in Θ using a numerical optimization algorithm to maximize the optimization objective function defined in step S2.

5. Step B4: Block-wise predict rigid parameters set {θ1, θ2,..., θ M} using sliding window strategy and fuse the rigid parameters set to output global rigid transformation matrix M final .

8. The medical image registration apparatus based on self-supervised representation learning and hierarchical feature optimization, characterized in that, Comprise: A pre-training module is configured to pre-train a 3D-PVT encoder using a mask autoencoder strategy to learn anatomical structure feature representation from unlabeled three-dimensional benchmark images. A registration module is configured to receive a benchmark image and a to-be-registered image of a patient, extract a multi-level feature pyramid through the pre-trained 3D-PVT encoder, and perform iterative optimization from coarse to fine layer by layer on the feature pyramid based on a mutual information maximization target, and output a registration result.

9. A computing device, characterized by Comprise: One or more processors; Memory; And one or more programs, wherein the one or more programs are stored in the memory and configured to be executed by one or more processors, one or more programs comprising instructions of the medical image registration method based on self-supervised representation learning and hierarchical feature optimization according to any one of the preceding claims 1-7.

10. Storage medium, characterized in that The storage medium stores one or more computer-readable programs, and the one or more programs comprise instructions suitable for being loaded by the memory and executing the medical image registration method based on self-supervised representation learning and hierarchical feature optimization according to any one of the preceding claims 1-7.

Images