A manifold learning-based multi-modal lower limb medical image registration method and system

Abstract

This application relates to the technical field of medical image processing, and more particularly to a multimodal lower limb medical image registration method and system based on manifold learning. The method first performs left-right separation on a whole-body lower limb CT scan to obtain independent left and right candidate regions; then extracts limb masks from local MRI and determines the side orientation. Local images are resampled using the overall resolution as a reference to achieve spatial resolution consistency. High-dimensional features of local blocks are extracted separately, a nearest neighbor graph is constructed, and the graph Laplacian characteristic equation is solved to obtain a multi-channel manifold feature map. Representational ambiguities in manifold embedding are eliminated through sign correction and orthogonal alignment, and then Fourier-Merlin transform is used to achieve coarse registration, followed by local refinement within a limited neighborhood. This invention completely separates and independently registers the left and right sides, avoiding interference from contralateral structures, and achieving stable localization and precise alignment of local MRI images in whole-body lower limb CT.

Description

Technical Field

[0001] This application relates to the technical field of medical image processing, and in particular to a multimodal lower limb medical image registration method and system based on manifold learning. Background Technology

[0002] In lower extremity disease assessment, preoperative planning, postoperative follow-up, and multimodal image fusion analysis, it is often necessary to align CT images covering the entire lower extremity with MRI images covering only the knee joint or adjacent regions of interest. CT has advantages in displaying bony structures, while MRI (Magnetic Resonance Imaging) has advantages in displaying soft tissues, bone marrow edema, and focal lesions; therefore, the combined use of the two has clear clinical value.

[0003] However, the registration between lower extremity local MRI and whole-body CT is not a traditional "complete image to complete image" registration, but rather a typical partial registration problem. The technical challenges mainly lie in the following three aspects: 1. Inconsistent anatomical coverage. Global CT usually includes extensive anatomical information of both lower limbs, while local MRI only covers one side or a local joint area. There is an inherent asymmetry in spatial coverage between the two, and traditional global registration methods cannot directly establish a stable correspondence.

[0004] 2. Spatial resolution is highly heterogeneous. Different modal images typically exhibit significant differences in voxel size, slice thickness, and anisotropy.

[0005] 3. Significant differences in body position. The lower limbs are prone to flexion, extension, external rotation, internal rotation, and relative displacement changes at different acquisition times, making it more difficult to locate local MRI in global CT. Local areas often cannot be accurately recovered through simple translation or single global optimization.

[0006] In existing technologies, for example, CN121482436A discloses a CT-MRI multimodal fusion intelligent grading method. Its core steps include: first, acquiring and standardizing CT and MRI images of the patient's brain; second, spatially aligning the two modalities using a multi-scale registration method; then, fusing the registered data and training an intelligent grading model using this fused data; and finally, automatically grading the glioma of a new patient using this model. However, this registration strategy relies on the assumption that the CT and MRI images have completely identical coverage areas. Traditional multi-scale registration methods cannot handle scenarios with severe mismatches in the field of view, such as "global-local" (e.g., whole lower limb CT and single knee MRI), and are prone to failure under interference from non-overlapping areas. Furthermore, this patent performs registration directly in the original or shallow feature space, lacking the ability to model the intrinsic structure of high-dimensional heterogeneous data, making it difficult to effectively bridge the significant modal differences between CT and MRI. For example, CN121221041A discloses an image fusion method that acquires three types of endoscopic images—white light, fluorescence, and infrared—and adaptively adjusts the focus; it uses different filtering methods for noise reduction, achieves spatial alignment based on feature point matching, and normalizes pixel values; it registers the registered image with preoperative CT or MRI images, then performs multi-scale decomposition, and reconstructs the fused image by adaptive weighted fusion; it extracts features such as color, texture, shape, and temperature, inputs them into a hybrid model to identify lesions and assess risk levels. However, this scheme registers white light, fluorescence, and infrared images acquired under the same endoscope, which is a homologous multimodal registration and cannot solve the problems of inconsistent coverage and lack of common anatomical areas between overall CT and local MRI; it uses a traditional feature point matching algorithm, which relies on pixel-level grayscale comparability and cannot handle registration scenarios where CT and MRI resolutions are heterogeneous and the physical meanings of grayscale differ. Summary of the Invention

[0007] The purpose of this invention is to provide a multimodal lower limb medical image registration method and system based on manifold learning, which partially solves or alleviates the registration difficulties between overall lower limb CT and local MRI caused by inconsistent anatomical coverage, highly heterogeneous spatial resolution, and large differences in body position in the prior art, and achieves stable positioning and accurate alignment of local images in the overall image.

[0008] To solve the aforementioned technical problems, the present invention specifically adopts the following technical solution: A first aspect of the present invention is to provide a multimodal lower limb medical image registration method based on manifold learning, comprising the following steps: S1: Acquire the first and second modal medical images to be registered. The first modal medical image is a CT image covering the entire range of both lower limbs, and the second modal medical image is an MRI image covering a local area. S2: Extract local limb masks from the second modality medical image and determine its side orientation; perform left-right separation on the first modality medical image and obtain the corresponding side candidate region corresponding to the side orientation of the second modality medical image; S3: Using the spatial resolution of the first modality medical image as a reference, resample the second modality medical image, and perform spatial resolution uniformity and preprocessing; S4: Extract high-dimensional feature vectors of local 3D neighborhood blocks centered on voxels from the corresponding candidate regions and the regions corresponding to the local limb masks, respectively, and construct a high-dimensional local block feature point set; construct a symmetric nearest neighbor graph based on the high-dimensional local block feature point set and determine the edge weights to form a weighted adjacency graph; construct a graph Laplacian matrix based on the weighted adjacency graph, solve for the low-dimensional embedding coordinates, and map them back to the original spatial position to obtain the first multi-channel manifold feature map and the second multi-channel manifold feature map; S5: Perform manifold alignment on the first multi-channel manifold feature map and the second multi-channel manifold feature map to eliminate the representation ambiguity caused by independent embedding of different modalities and obtain a set of aligned manifold feature maps; S6: Input the aligned set of manifold feature maps into the frequency domain medical image registration module, and estimate the scale factor, rotation parameter and translation parameter based on Fourier-Merlin transform to obtain the coarse medical image registration transformation; S7: Based on the coarse medical image registration transformation, a restricted neighborhood search and refinement optimization are performed on the overlapping area between the region corresponding to the local limb mask and the corresponding side candidate region to obtain the final spatial transformation. S8: Apply the final spatial transformation to the second modality medical image, map it to the spatial coordinate system of the first modality medical image, and output the registration result.

[0009] Furthermore, in step S2, performing left-right separation on the first modality medical image includes: Extract the lower limb surface mask from the first modality of medical images; If the lower limb surface mask contains two main connected regions, then the left and right are marked according to the horizontal coordinate of the physical centroid. The candidate region with the larger horizontal coordinate is the left candidate region, and the candidate region with the smaller horizontal coordinate is the right candidate region. If the lower limb surface mask is a single connected region, then perform bi-cluster clustering on the two-dimensional projection point set of the lower limb surface mask, estimate the separating polylines in the background gaps between the bi-clusters, and perform separation based on the separating polylines.

[0010] Furthermore, in step S2, extracting a local limb mask from the second modality medical image and determining its side includes: Extract local limb masks from the second modality of medical images; Based on the horizontal coordinates of the physical centroid of the connected domain of the local limb mask or the morphological matching score of the local limb mask, the second modality medical image is determined to be left, right, or bilateral. The morphological matching score is calculated as a weighted sum of the following indicators based on the local limb mask: The similarity of shapes, consistency of physical dimensions, consistency of volume ratios, and consistency of center position after alignment in the main direction.

[0011] Furthermore, in step S3, the resampling adopts three-dimensional linear interpolation, constructs a resampling grid based on the voxel spacing of the first modality medical image, and constrains the target voxel size after resampling based on the principle of maintaining physical range consistency.

[0012] Furthermore, in step S4, the edge weights of the nearest neighbor graph are defined as follows: When two sampling points are neighbors, the edge weight is equal to the value of an exponential function with the square of the difference between the eigenvectors of the two points divided by the negative of the product of the local scale parameters of the two points; otherwise, the edge weight is zero. The local scale parameter of each sampling point is defined as the average distance between the current point and all its neighbors.

[0013] Furthermore, in step S5, performing manifold alignment on the first multi-channel manifold feature map and the second multi-channel manifold feature map includes: Let the embedding matrix of the corresponding candidate region be the first embedding matrix, and the embedding matrix of the region corresponding to the local limb mask be the second embedding matrix. First, determine the symbol correction vector based on the correlation between each channel. Each component of the symbol correction vector is positive or negative. Then, solve for an orthogonal alignment matrix such that the Frobenius norm of the difference between the first embedding matrix and the second embedding matrix multiplied by the symbol correction diagonal matrix and then multiplied by the orthogonal alignment matrix is ​​minimized, thus obtaining a set of aligned manifold feature maps.

[0014] Furthermore, the aligned set of manifold feature maps includes a first multi-channel manifold feature map that retains the original manifold embedding, and an aligned manifold feature map after performing a manifold alignment transformation on the second multi-channel manifold feature map; The aligned manifold feature map is obtained by multiplying the second embedding matrix by the sign-corrected diagonal matrix and the orthogonal alignment matrix on the right in sequence.

[0015] Furthermore, in step S6, based on the Fourier-Merlin transform, the scale factor, rotation parameters, and translation parameters are estimated to obtain the coarse medical image registration transform, including: The aligned set of manifold feature maps are unified into a common computational grid. The scaling factor is estimated based on the radial distribution of the spectrum, the rotation parameters are estimated based on the multi-level angle search of the spectrum amplitude, and the translation parameters are determined by combining phase correlation and translation candidate evaluation.

[0016] Furthermore, S7 includes: Obtain the body surface mask corresponding to the candidate region on the corresponding side and construct a distance map. Optimize within the preset translation and local pose perturbation range near the coarse medical image registration transformation. Only when the optimization result increases the overlap between the local limb mask and the body surface mask corresponding to the candidate region on the corresponding side or reduces the center offset, and the optimization displacement does not exceed the local drift threshold, is the optimization result accepted as the final spatial transformation.

[0017] Secondly, this application also discloses a multimodal lower limb medical image registration system based on manifold learning, the system comprising: The medical image acquisition module is configured to acquire a first modality medical image and a second modality medical image to be registered. The first modality medical image is a CT image covering the entire range of both lower limbs, and the second modality medical image is an MRI image covering a local area. The lateral identification module is configured to extract local limb masks from the second modality medical image and determine their laterality, perform left-right separation on the first modality medical image, and obtain the corresponding side candidate region corresponding to the laterality of the second modality medical image; The resampling module is configured to resample the second modality medical image with reference to the spatial resolution of the first modality medical image, and to perform spatial resolution uniformity and preprocessing. The manifold feature map acquisition module is configured to extract high-dimensional feature vectors of local three-dimensional neighborhood blocks centered on voxels from the corresponding candidate regions on the sides and the regions corresponding to the local limb masks, respectively, and construct a high-dimensional local block feature point set; construct a symmetric nearest neighbor graph based on the high-dimensional local block feature point set and determine the edge weights to form a weighted adjacency graph; construct a graph Laplacian matrix based on the weighted adjacency graph, solve for the low-dimensional embedding coordinates, and map them back to the original spatial position to obtain a first multi-channel manifold feature map and a second multi-channel manifold feature map; The manifold alignment module is configured to perform manifold alignment on the first multi-channel manifold feature map and the second multi-channel manifold feature map to eliminate representational ambiguity caused by independent embedding of different modalities and obtain a set of aligned manifold feature maps. The coarse registration transformation module is configured to input the aligned set of manifold feature maps into the frequency domain medical image registration module, and estimate the scale factor, rotation parameters and translation parameters based on the Fourier-Merlin transform to obtain the coarse medical image registration transformation; The fine registration transformation module is configured to perform restricted neighborhood search and refinement optimization on the overlapping area of ​​the region corresponding to the local limb mask and the corresponding side candidate region, based on the coarse medical image registration transformation, to obtain the final spatial transformation. The registration module is configured to apply the final spatial transformation to the second modality medical image, map it to the spatial coordinate system of the first modality medical image, and output the registration result.

[0018] Beneficial technical effects: 1. This invention addresses the inconsistency between the anatomical coverage of whole-limb CT and local MRI by employing a strategy of complete separation and independent registration of the left and right sides. By separating the whole-limb CT scan to obtain independent left and right candidate regions, and simultaneously extracting a limb mask from the local MRI and determining the side, the local MRI is registered only with the ipsilateral candidate region of the whole-limb CT. This strategy fundamentally avoids interference from contralateral limb structures, preventing local images from being incorrectly mapped to the contralateral limb during global registration, and achieving stable ipsilateral localization of the local MRI within the whole-limb CT. Even if the local MRI only covers a local area such as the knee, hip, or ankle joint on one side, this invention can still accurately register it to the corresponding anatomical location in the whole-limb CT.

[0019] 2. This invention addresses the issue of flexion-extension, external rotation, internal rotation, and relative displacement changes in the lower limbs at different acquisition times by designing a multi-level registration strategy. First, coarse registration is achieved in the frequency domain using Fourier-Merlin transform to estimate the scale factor, rotation parameters, and translation parameters. Then, local refinement optimization is performed within a restricted neighborhood near the coarse registration transform. By constructing a surface mask distance map and optimizing mask overlap and center offset, registration accuracy is further improved. This multi-level strategy ensures both robust registration and sub-voxel-level precise alignment.

[0020] 3. To address the representational ambiguity caused by independent manifold embedding for different modes, this invention designs an alignment method combining sign correction and orthogonal alignment. First, the sign correction vector is determined based on the correlation between channels to eliminate ambiguity in sign direction. Then, the orthogonal alignment matrix is ​​solved to eliminate rotational or reflection differences between the low-dimensional embedding bases of different modes. Through this alignment operation, the manifold feature maps of the two modes achieve unified representation in the same low-dimensional feature space, providing comparable feature representations for subsequent frequency domain registration.

[0021] 4. This invention addresses the problem of large differences in voxel spacing between CT and MRI. It resamples local MRI data using the overall spatial resolution of CT as a reference, ensuring consistent or nearly consistent voxel spacing between the two modalities in physical space. Based on this, high-dimensional features of local 3D neighborhood blocks are extracted and dimensionality is reduced using manifold learning methods to obtain a multi-channel manifold feature map. This manifold feature map effectively captures structural information within the local neighborhood while reducing the dependence on direct comparability of original grayscale values, thus eliminating the impact of spatial resolution differences on feature similarity measurement. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. In all the drawings, similar elements or parts are generally identified by similar reference numerals. The elements or parts in the drawings are not necessarily drawn to scale. Obviously, the drawings described below are some embodiments of the present invention, and those skilled in the art can obtain other drawings based on these drawings without any creative effort.

[0023] Figure 1 This is a flowchart of a multimodal lower limb medical image registration method based on manifold learning in one embodiment of this application.

[0024] Figure 2 This is a schematic diagram of left-right separation of the overall lower limb image and side separation determination of the local image in one embodiment of this application.

[0025] Figure 3a The manifold feature map of channel 1 for the overall candidate region of CT.

[0026] Figure 3b The manifold feature map of channel 2 for the overall candidate region of CT.

[0027] Figure 3c The manifold feature map of channel 3 for the overall candidate region of CT.

[0028] Figure 3d Manifold feature map of channel 1 of the aligned MRI local image.

[0029] Figure 3e Manifold feature map of channel 2 in the aligned MRI local image.

[0030] Figure 3f Manifold feature map of channel 3 in the aligned local MRI image.

[0031] Figure 4a A schematic diagram of the registration results sliced ​​along the axial direction.

[0032] Figure 4b A schematic diagram of the coronal section showing the registration results.

[0033] Figure 4c A schematic diagram of the registration results in the sagittal plane.

[0034] Figure 5 This is a schematic diagram of the module structure of a multimodal lower limb medical image registration system based on manifold learning in one embodiment of this application. Detailed Implementation

[0035] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0036] In this document, suffixes such as "module," "part," or "unit" used to denote elements are used only for the purpose of illustrative purposes and have no specific meaning in themselves. Therefore, "module," "part," or "unit" may be used interchangeably.

[0037] In this document, the terms "upper," "lower," "inner," "outer," "front," "rear," "one end," and "the other end," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing the present invention and for simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on the present invention. Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0038] In this document, unless otherwise explicitly specified and limited, the terms "installed," "equipped with," "connected," etc., should be interpreted broadly. For example, "connection" can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection, a direct connection, or an indirect connection through an intermediate medium; it can be a connection within two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0039] In this document, "and / or" includes any and all combinations of one or more of the listed related items.

[0040] In this article, "multiple" means two or more, that is, it includes two, three, four, five, etc.

[0041] In this paper, "highly heterogeneous spatial resolution" is used to indicate that the two modes differ significantly in terms of sampling density and layer resolution. This difference directly affects the comparability of local structural features and the stability of similarity measures.

[0042] Figure 1 This paper presents a flowchart of a multimodal lower limb medical image registration method based on manifold learning, as described in this application. (Refer to...) Figure 1 The method specifically includes the following steps: S1: Acquire a first modal medical image and a second modal medical image to be registered. The first modal medical image is a CT image covering the entire range of both lower limbs, and the second modal medical image is an MRI image covering a local area.

[0043] In some embodiments, the first modality of medical image is a CT image covering the entire extent of both lower limbs. For example, the CT image extends from the upper edge of the pelvis to below the ankle, fully encompassing the skeletal, joint, and soft tissue structures of the left and right lower limbs, including the hip, femur, knee, tibia, fibula, and ankle joints. This CT image is used to provide a global anatomical reference space as a fixed reference image for subsequent medical image registration (also known as partial registration).

[0044] The second modality of medical image is an MRI image covering a localized region. For example, this MRI image is an MRI image covering a region of interest in a single lower limb, including but not limited to the knee, hip, or ankle joints. This MRI image is used to provide localized high-resolution soft tissue information as a floating image for subsequent medical image registration.

[0045] It should be noted that the first and second modal medical images are inconsistent in terms of anatomical coverage. The first modal medical images cover both complete lower limbs, while the second modal medical images only cover a local area on one side. They are also heterogeneous in terms of spatial resolution. The voxel spacing of the first modal medical images is usually smaller, while the voxel spacing of the second modal medical images is usually larger. There may also be differences in the acquisition position between the two, including but not limited to changes in the degree of flexion and extension, internal and external rotation angles, and relative displacement of the lower limbs.

[0046] S2: Extract local limb masks from the second modality medical image and determine its side orientation; perform left-right separation on the first modality medical image and obtain the corresponding side candidate region corresponding to the side orientation of the second modality medical image.

[0047] In some embodiments, step S2, extracting a local limb mask from the second modality medical image and determining its side, includes: Extract local limb masks from the second modality of medical images; Based on the horizontal coordinates of the physical centroid of the connected domain of the local limb mask or the morphological matching score of the local limb mask, the second modality medical image is determined to be left, right, or bilateral. The morphological matching score is calculated as a weighted sum of the following indicators based on the local limb mask: The similarity of shapes, consistency of physical dimensions, consistency of volume ratios, and consistency of center position after alignment in the main direction.

[0048] In one specific implementation, threshold segmentation is performed on the second modality medical image to separate the human tissue from the background air, obtaining an initial binary mask. Then, connected component analysis is performed on the initial binary mask to select the connected component with the largest area as the main limb region, resulting in a local limb mask. This local limb mask is used to characterize the area where the limb is located in the second modality medical image.

[0049] Then, based on the horizontal coordinates of the physical centroid of the connected components of the local limb mask, or based on the morphological matching score of the local limb mask, it is determined whether the second modality medical image belongs to the left, right, or bilateral side. Specifically, the horizontal coordinate value of the physical centroid of the connected components of the local limb mask is calculated. If the horizontal coordinate value is significantly biased to one side, the side is directly determined based on the horizontal coordinate of the physical centroid. If the horizontal coordinate of the physical centroid of the local limb mask is located in the middle region, the morphological matching score of the local limb mask is further calculated. This morphological matching score is compared with the left standard template and the right standard template respectively. The side corresponding to the higher matching score is selected as the determination result. If the matching score with the bilateral template (which includes templates from both sides) exceeds a preset threshold, it is determined to be bilateral. If the morphological matching score of the local limb mask is equal to or similar to the left and right standard templates, and the matching score with both templates does not exceed the threshold, it is determined to be an invalid second modality medical image. (When the morphological matching score of the local limb mask is equal to or similar to the left and right standard templates, it indicates that the mask has achieved a considerable degree of matching with both templates, meaning that the mask simultaneously includes the morphology of both sides. Under normal circumstances, the second modality medical image can be determined to belong to the bilateral category (i.e., a high degree of matching with both templates). Therefore, under normal circumstances, it is unlikely that the morphological matching score of the local limb mask is equal to or similar to the left and right standard templates, and the matching score with both templates does not exceed the preset threshold. If this does occur, it indicates that the image is an invalid second modality medical image.) Alternatively, if the morphological matching score of the local limb mask is zero when compared with both the left and right standard templates, it also indicates that the image is an invalid second modality medical image. In this case, an invalid second modality medical image should be used instead, or manual intervention should be performed for determination.

[0050] The morphological matching score is calculated as a weighted sum of the following indicators for the local limb mask: shape similarity after main direction alignment, physical size consistency, volume ratio consistency, and center position consistency. Specifically, shape similarity after main direction alignment refers to the shape similarity calculated after rotating the local limb mask to align with the main direction of the standard template; physical size consistency refers to the degree of matching between the physical span of the local limb mask and the standard template in each direction; volume ratio consistency refers to whether the ratio of the overall volume of the local limb mask to the standard template meets expectations; and center position consistency refers to the spatial distance between the physical centroid of the local limb mask and the centroid of the standard template. The morphological matching score is obtained by assigning preset weight coefficients to each of the above four indicators and summing them. For example, the weight of shape similarity is 0.4, and the weights of physical size consistency, volume ratio consistency, and center position consistency are each 0.2, meaning shape similarity is given a higher weight, while the other three indicators have equal weights. The calculated morphological matching score is compared with the left and right standard templates respectively, and the side corresponding to the higher score is taken as the judgment result.

[0051] In some embodiments, performing left-right separation on the first modal medical image in step S2 includes: Extract the lower limb surface mask from the first modality of medical images; If the lower limb surface mask contains two main connected regions, then the left and right are marked according to the horizontal coordinate of the physical centroid. The candidate region with the larger horizontal coordinate is the left candidate region, and the candidate region with the smaller horizontal coordinate is the right candidate region. If the lower limb surface mask is a single connected region, then perform bi-cluster clustering on the two-dimensional projection point set of the lower limb surface mask, estimate the separating polylines in the background gaps between the bi-clusters, and perform separation based on the separating polylines.

[0052] In one specific implementation, a lower limb surface mask is extracted from the first modality of medical images. Specifically, by performing body surface thresholding on the first modality of medical images, the human tissue in the image is separated from the background air to obtain an initial binary mask; then, connected component analysis is performed on the initial binary mask to select one or more connected components with the largest area as the lower limb surface mask. The lower limb surface mask is used to characterize the body surface contour range of both lower limbs in the first modality of medical images.

[0053] If the lower limb surface mask contains two main connected regions, the physical centroid of each connected region is calculated, and the lateral coordinates of the two centroids are compared. The region with the larger lateral coordinate is designated as the left candidate region, and the region with the smaller lateral coordinate is designated as the right candidate region. Specifically, let the two main connected regions be C1 and C2, and their physical centroids be μ(C1) and μ(C2), respectively. The lateral coordinate values ​​of μ(C1) and μ(C2) are compared, and the connected region with the larger lateral coordinate value is designated as the left candidate region, and the connected region with the smaller lateral coordinate value is designated as the right candidate region.

[0054] If the lower limb surface mask is a single connected region, meaning that both lower limbs are connected as one unit in the surface mask and cannot be directly separated, then bi-cluster clustering is performed on the two-dimensional projection point set of the lower limb surface mask. A separating polyline is estimated in the background gap between the two clusters, and separation is performed based on the separating polyline. Specifically, the lower limb surface mask is projected onto a horizontal plane to obtain a two-dimensional projection point set; bi-cluster clustering is performed on the two-dimensional projection point set to divide it into two cluster subsets; in the background gap between the two cluster subsets, a separating polyline is estimated by optimizing the background cost function and smoothing constraint term, and then the lower limb surface mask is segmented into a left candidate region and a right candidate region based on the separating polyline. Using the above method, even if the lower limbs appear adhered in the body surface mask (i.e., the lower limb mask is a single connected region; for example, in the extracted binary image of the lower limb mask, the regions of the left and right legs are connected, appearing adhered, and the adhered part in the middle does not belong to any cluster after clustering, making it difficult to distinguish directly), effective separation can still be achieved. Specifically, after clustering determines the approximate left and right regions, by optimizing the background cost function and smoothing constraint term, a continuous and smooth dividing line is estimated between the two cluster subsets. Based on this dividing line, the adhered body surface mask is cut into independent left and right candidate regions, thereby achieving effective separation in the adhered state.

[0055] In some embodiments, a background dividing polyline is defined. The optimization formula is: ; in, Represents the background cost function. The smoothing constraint term is represented by minimizing the above formula to obtain a dividing polyline that passes through the real background gap between the two lower limbs as much as possible.

[0056] In some embodiments, the background cost function may take the following form: ; in, The dividing line is located on the vertical axis. The x-coordinate value at that location, These are the candidate points that the broken line passes through; and These represent the candidate point to the left cluster subset, respectively. and right-side cluster subset The shortest Euclidean distance, Take the smaller of the two; This is a normalization constant, which can be the distance between the centroids of two cluster subsets or a fixed parameter. This cost function is only applicable when... and The optimization is performed in the middle region. When the candidate point is located in the middle of the two clusters, the background cost function is minimized, thus guiding the separating line to preferentially pass through the middle of the two clusters.

[0057] In the smoothing constraint term, λ is the smoothing constraint coefficient, and in this embodiment, the smoothing constraint coefficient λ is set to 2.0. The smoothing constraint term penalizes the change in the x-coordinate of the polyline between adjacent rows, forcing the separating polyline to maintain a smooth transition and avoiding violent swings or jagged shapes, thereby ensuring that the separation boundary conforms to the smoothness expectation of the lower limb anatomy.

[0058] S3: Using the spatial resolution of the first modality medical image as a reference, resample the second modality medical image, and perform spatial resolution uniformity and preprocessing.

[0059] In step S3, the resampling uses three-dimensional linear interpolation to construct a resampling grid based on the voxel spacing of the first modality medical image, and constrains the size of the target voxel after resampling based on the principle of maintaining physical range consistency.

[0060] In one specific implementation, the voxel spacing of the first modality medical image is assumed to have three values ​​in the left-right, front-back, and head-to-toe directions, while the original voxel spacing of the second modality medical image has three other values ​​in the same three directions. Due to differences in acquisition parameters between the two modalities, their voxel spacings are usually not equal; for example, CT images have smaller voxel spacings while MRI images have larger voxel spacings. This high heterogeneity in spatial resolution directly affects the stability of subsequent feature extraction and similarity measurement. To eliminate this difference, the second modality medical image is resampled using the voxel spacing of the first modality medical image as a reference.

[0061] In some embodiments, resampling employs three-dimensional linear interpolation, constructing a resampling grid based on the voxel spacing of the first modality medical image, and constraining the target voxel size after resampling based on the principle of maintaining physical range consistency. Let the original voxel size of the second modality medical image be the number of voxels contained in the left-right, front-back, and head-to-toe directions, respectively. The target voxel size after resampling is determined according to the principle of maintaining physical range consistency, i.e., the target voxel size in each direction is equal to the original voxel size multiplied by the ratio of the second modality voxel spacing to the first modality voxel spacing in that direction. Subsequently, a resampling grid is constructed using the voxel spacing of the first modality medical image as the target voxel spacing. For each target voxel in the resampling grid, its coordinates are located in its corresponding original local image space, and the grayscale value of the target voxel is calculated using three-dimensional linear interpolation. For target voxels outside the original local image domain, their grayscale values ​​are assigned a preset background value. Through the above resampling operation, a second modality medical image with the same spatial resolution as the first modality medical image is obtained.

[0062] Preferably, after resampling, preprocessing operations can be performed on the resampled second modality medical image. The preprocessing includes, but is not limited to: smoothing the image to reduce noise, normalizing the image intensity to eliminate grayscale differences caused by different acquisition parameters, and suppressing artifact regions in the image. After resampling and preprocessing, the second modality medical image achieves the same spatial resolution as the first modality medical image, laying the foundation for subsequent manifold feature extraction and registration operations.

[0063] S4: Extract high-dimensional feature vectors of local 3D neighborhood blocks centered on voxels from the corresponding candidate regions on the sides and the regions corresponding to the local limb masks, respectively, and construct a high-dimensional local block feature point set; construct a symmetric nearest neighbor graph based on the high-dimensional local block feature point set and determine the edge weights to form a weighted adjacency graph; construct a graph Laplacian matrix based on the weighted adjacency graph, solve for the low-dimensional embedding coordinates, and map them back to the original spatial position to obtain the first multi-channel manifold feature map and the second multi-channel manifold feature map.

[0064] In step S4, the edge weights of the nearest neighbor graph are defined as follows: When two sampling points are neighbors, the edge weight is equal to the value of an exponential function with the square of the difference between the eigenvectors of the two points divided by the negative of the product of the local scale parameters of the two points; otherwise, the edge weight is zero. The local scale parameter of each sampling point is defined as the average distance between the current point and all its neighbors.

[0065] In one specific implementation, for each voxel in the corresponding candidate region, a local 3D neighborhood block is extracted centered on that voxel, and the gray values ​​of all voxels within that neighborhood block are expanded into a high-dimensional feature vector. Similarly, for each voxel in the region corresponding to the local limb mask, a local 3D neighborhood block is also extracted centered on that voxel and expanded into a high-dimensional feature vector. The high-dimensional feature vectors corresponding to all voxels together constitute a high-dimensional local block feature point set.

[0066] In some embodiments, the edge weight of the nearest neighbor graph is defined as follows: when two sampled points are neighbors, the edge weight is equal to the value of an exponential function with the square of the Euclidean distance between the difference of the feature vectors of the two points divided by the negative of the product of the local scale parameters of the two points; otherwise, the edge weight is zero. The local scale parameter of each sampled point is defined as the average distance between the current point and all its neighbors. Specifically, for any two sampled points, it is first determined whether they are neighbors, i.e., whether each point is in the neighbor set of the other. If the condition of being neighbors is met, the Euclidean distance between the feature vectors of the two sampled points is calculated, and the square of this distance is taken. Then, the local scale parameters of the two sampled points are calculated separately, and the local scale parameter of each sampled point is defined as the average Euclidean distance between that point and all its neighbors. The square of the Euclidean distance is divided by the product of the two local scale parameters, and the result is the edge weight. If the condition of being neighbors is not met, the edge weight is directly set to zero. By traversing all pairs of sampled points in the above manner, a weighted adjacency matrix is ​​obtained.

[0067] Based on the weighted adjacency matrix, a corresponding degree matrix is ​​constructed, where the elements on the diagonal of the degree matrix are the sum of the edge weights of all edges in the corresponding row of the weighted adjacency matrix, and the elements outside the diagonal are zero. Then, a graph Laplacian matrix is ​​constructed, which is equal to the degree matrix minus the weighted adjacency matrix. The generalized characteristic equation is solved to obtain eigenvectors corresponding to several non-trivial eigenvalues; these eigenvectors are the low-dimensional embedding coordinates of the high-dimensional local block feature point set. The low-dimensional embedding coordinate components of each sampling point are mapped back to the corresponding voxels according to their positions in the original voxel space, forming a multi-channel manifold feature map. The first multi-channel manifold feature map is obtained by extracting and processing the corresponding candidate region, and the second multi-channel manifold feature map is obtained by extracting and processing the region corresponding to the local limb mask.

[0068] S5: Perform manifold alignment on the first multi-channel manifold feature map and the second multi-channel manifold feature map to eliminate representational ambiguity caused by independent embedding of different modalities and obtain a set of aligned manifold feature maps.

[0069] In step S5, performing manifold alignment on the first multi-channel manifold feature map and the second multi-channel manifold feature map includes: Let the embedding matrix of the corresponding candidate region be the first embedding matrix, and the embedding matrix of the region corresponding to the local limb mask be the second embedding matrix. First, determine the symbol correction vector based on the correlation between each channel. Each component of the symbol correction vector is positive or negative. Then, solve for an orthogonal alignment matrix such that the Frobenius norm of the difference between the first embedding matrix and the second embedding matrix multiplied by the symbol correction diagonal matrix and then multiplied by the orthogonal alignment matrix is ​​minimized, thus obtaining a set of aligned manifold feature maps.

[0070] In some embodiments, the aligned set of manifold feature maps includes a first multi-channel manifold feature map that retains the original manifold embedding, and an aligned manifold feature map after performing a manifold alignment transformation on the second multi-channel manifold feature map.

[0071] The aligned manifold feature map is obtained by multiplying the second embedding matrix by the sign-corrected diagonal matrix and the orthogonal alignment matrix on the right in sequence.

[0072] In some embodiments, since the first modality medical image and the second modality medical image are embedded in manifolds independently, there is a problem of non-uniqueness in their representations in the low-dimensional feature space. Specifically, on the one hand, the feature vector of each embedded channel can be multiplied by negative one without changing the manifold structure, resulting in sign ambiguity; on the other hand, there may be differences in orthogonal bases such as rotation or reflection between the low-dimensional embedding bases obtained from different modalities, causing inconsistent directional representations of the same anatomical structure in the low-dimensional coordinate system. To solve the above problems, step S5 performs a manifold alignment operation on the two manifold feature maps.

[0073] In step S5, performing manifold alignment on the first and second multi-channel manifold feature maps includes: setting the embedding matrix of the corresponding candidate region as the first embedding matrix, and the embedding matrix of the region corresponding to the local limb mask as the second embedding matrix. The number of rows in the first embedding matrix corresponds to the number of sampling points in the corresponding candidate region, and the number of columns corresponds to the number of manifold embedding channels; the number of rows in the second embedding matrix corresponds to the number of sampling points in the region corresponding to the local limb mask, and the number of columns also corresponds to the number of manifold embedding channels.

[0074] First, a symbol correction vector is determined based on the correlation between the channels of the first and second embedding matrices. Each component of this symbol correction vector takes a value of either positive or negative one, and is used to correct the symbol direction of each channel, so that the embedding directions of the two modes in the corresponding channels tend to be consistent. Specifically, for each channel, the correlation between the column vector of that channel in the first embedding matrix and the column vector of that channel in the second embedding matrix is ​​calculated. If the correlation is positive, the symbol correction component takes a value of positive one; if the correlation is negative, it takes a value of negative one; if the correlation is close to zero, it can be randomly assigned a value of positive or negative one.

[0075] Secondly, based on the sign correction, an orthogonal alignment matrix is ​​obtained. This orthogonal alignment matrix satisfies the constraint that its transpose multiplied by itself equals the identity matrix; that is, the matrix is ​​orthogonal. The optimization objective is to minimize the Frobenius norm of the difference between the first and second embedding matrices, which are then right-multiplied successively by the sign-corrected diagonal matrix and then right-multiplied by the orthogonal alignment matrix. The Frobenius norm measures the overall magnitude of the difference between two matrices; by minimizing this norm, the transformed second embedding matrix is ​​made as close as possible to the first embedding matrix. The solution to this optimization problem can be obtained through orthogonal Protodyakonov analysis.

[0076] Through the aforementioned sign correction and orthogonal alignment steps, a set of aligned manifold feature maps is obtained. This set includes a first multi-channel manifold feature map that retains the original manifold embedding, and an aligned manifold feature map obtained by performing a manifold alignment transformation on the second multi-channel manifold feature map. The first multi-channel manifold feature map remains unchanged throughout the alignment process, serving as a fixed reference. The second multi-channel manifold feature map, after sign correction and orthogonal alignment transformation, is mapped to the low-dimensional feature space containing the first multi-channel manifold feature map. The aligned manifold feature map is obtained by multiplying the second embedding matrix sequentially by the sign correction diagonal matrix and the orthogonal alignment matrix, and then mapping this result back to the corresponding voxel positions according to the original voxel space positions. Thus, the manifold feature maps of the two modes achieve unified representation in the same low-dimensional feature space, providing comparable feature representations for subsequent frequency domain registration.

[0077] S6: Input the aligned set of manifold feature maps into the frequency domain medical image registration module, and estimate the scale factor, rotation parameter and translation parameter based on Fourier-Merlin transform to obtain the coarse medical image registration transformation.

[0078] In step S6, based on the Fourier-Merlin transform, the scale factor, rotation parameters, and translation parameters are estimated to obtain the coarse medical image registration transform, which includes: The aligned set of manifold feature maps are unified into a common computational grid. The scaling factor is estimated based on the radial distribution of the spectrum, the rotation parameters are estimated based on the multi-level angle search of the spectrum amplitude, and the translation parameters are determined by combining phase correlation and translation candidate evaluation.

[0079] In one specific implementation, after the manifold alignment in step S5, the manifold feature maps of the first modality medical image and the second modality medical image have achieved unified representation in the low-dimensional feature space, but there are still differences in scale, rotation, and translation between the two in physical space. In order to quickly obtain the initial spatial correspondence, step S6 uses a frequency domain method for coarse registration.

[0080] In step S6, the coarse medical image registration transformation is obtained by estimating the scale factor, rotation parameters, and translation parameters based on the Fourier-Merlin transform, including: unifying a set of aligned manifold feature maps into a common computation grid, estimating the scale factor based on the spectral radial distribution, estimating the rotation parameters based on the multi-level angle search of the spectral amplitude, and determining the translation parameters by combining phase correlation and translation candidate evaluation.

[0081] Specifically, the aligned first and second multi-channel manifold feature maps are first transformed onto the same computational grid. Since the two manifold feature maps may have different sizes or sampling intervals, they need to be interpolated into a unified grid for frequency domain analysis. After unifying the grid, Fourier transforms are performed on the two manifold feature maps respectively, transforming them from the spatial domain to the frequency domain.

[0082] In the frequency domain, scaling manifests as radial scaling of the spectral amplitude. By transforming the spectral amplitude to log-polar coordinates, the scale factor and rotation angle can be represented as translations in radial and angular coordinates, respectively. Specifically, the radial distribution of the spectral amplitude is calculated, and the scale factor between two images is estimated by performing a one-dimensional correlation analysis on the radial projection.

[0083] For the estimation of rotation parameters, a multi-level angle search strategy based on spectral amplitude is adopted. First, a global search is performed at coarser angle intervals to initially determine the candidate range of rotation angles; then, a local search is performed within the candidate range at finer angle intervals to obtain the precise rotation parameters. This multi-level search strategy can reduce computational complexity while ensuring estimation accuracy.

[0084] After determining the scale and rotation parameters, the original manifold feature map is scaled and rotated for correction, at which point the only difference between the two images is translation. The phase correlation method is used to calculate the translation parameters between the two images. This method utilizes the phase information of the cross-power spectrum in the frequency domain, enabling robust estimation of translation amounts at the integer pixel level. To further improve the accuracy of the translation parameters, multiple translation candidates can be generated near the integer translation estimation results, and the optimal translation parameter is selected by combining the correlation evaluation index in the spatial domain.

[0085] Through the frequency domain estimation described above, the scale factor, rotation matrix, and translation vector are obtained, which together constitute a coarse medical image registration transformation. This transformation can roughly align the second modality medical image to the spatial coordinate system of the first modality medical image, providing good initial values ​​for the local thinning in step S7.

[0086] S7: Based on the coarse medical image registration transformation, a restricted neighborhood search and refinement optimization are performed on the overlapping area between the region corresponding to the local limb mask and the corresponding side candidate region to obtain the final spatial transformation.

[0087] In some embodiments, S7 includes: Obtain the body surface mask corresponding to the candidate region on the corresponding side and construct a distance map. Optimize within the preset translation and local pose perturbation range near the coarse medical image registration transformation. Only when the optimization result increases the overlap between the local limb mask and the body surface mask corresponding to the candidate region on the corresponding side or reduces the center offset, and the optimization displacement does not exceed the local drift threshold, is the optimization result accepted as the final spatial transformation.

[0088] In some embodiments, after the frequency domain coarse registration in step S6, the second modality medical image has been initially aligned to the spatial coordinate system of the first modality medical image. However, due to the limitations of the frequency domain method and the nonlinear differences between the two modalities, the coarse registration result may have local deviations. Step S7 performs local search and optimization within a small range near the coarse registration transformation to further improve the registration accuracy.

[0089] Specifically, firstly, a surface mask for the corresponding candidate region is obtained from the first modality of medical image to represent the surface contour of the corresponding lower limb. A distance map is then constructed based on this surface mask. Each voxel position in the distance map stores the shortest distance from that position to the boundary of the surface mask; a positive distance value indicates the location is inside the mask, a negative distance value indicates the location is outside the mask, and a distance value of zero indicates the location is on the boundary. This distance map serves as a reference for subsequent optimization of the objective function.

[0090] Secondly, a restricted search space is defined for the vicinity of the coarse medical image registration transformation. This restricted search space includes a translational perturbation range and a local pose perturbation range. The translational perturbation range refers to the preset offset interval allowed in each of the left-right, front-back, and head-to-toe directions, based on the coarse registration translation parameters; for example, the offset in each direction shall not exceed a certain number of millimeters. The local pose perturbation range refers to the preset small angle rotation allowed around each coordinate axis based on the coarse registration rotation parameters; for example, the rotation in each direction shall not exceed a certain number of degrees. The setting of this restricted search space ensures sufficient degrees of freedom for refinement optimization while avoiding getting trapped in local extrema far from the true solution.

[0091] Then, an optimization search is performed within the restricted search space. For each candidate transformation, the local limb mask is transformed according to the candidate transformation, and the overlap between the transformed local limb mask and the corresponding side candidate region surface mask, as well as the spatial distance between their physical centroids, are calculated. The overlap can be measured by the Dice coefficient or the Jaccard coefficient; a higher overlap indicates better alignment, and a smaller center offset indicates closer positions.

[0092] During the optimization process, multiple sets of preset translation and local pose perturbation parameters are used to generate several candidate transformations. For each candidate transformation, its corresponding overlap and center offset are calculated and compared with the result of the coarse registration transformation. Only when a candidate transformation can increase the overlap or decrease the center offset, and the displacement of the candidate transformation relative to the coarse registration transformation does not exceed a preset local drift threshold, is the candidate transformation accepted as the optimization result.

[0093] Through the aforementioned restricted neighborhood search and refinement optimization, the final spatial transformation is obtained. This final transformation performs local fine-tuning adjustments based on the coarse registration result, ensuring improved registration accuracy while avoiding misalignment caused by excessive deviation from the coarse registration result through drift threshold constraints. Applying the final spatial transformation to the second modality of medical images enables precise localization of local MRI within the overall CT scan.

[0094] S8: Apply the final spatial transformation to the second modality medical image, map it to the spatial coordinate system of the first modality medical image, and output the registration result.

[0095] In some embodiments, the final spatial transformation obtained in step S7 includes a scale factor, a rotation matrix, and a translation vector. For each voxel in the second modality medical image, the physical coordinates of the voxel in the local image space are first obtained, and then the corresponding position of the voxel in the first modality medical image spatial coordinate system is calculated through the final spatial transformation. Since the transformed coordinates may fall on non-integer grid points, an appropriate interpolation method, such as trilinear interpolation or spline interpolation, is used to calculate the gray value at that position. All voxels in the second modality medical image are traversed to generate a transformed image located in the first modality medical image spatial coordinate system.

[0096] When outputting the registration results, several formats can be provided. One format is to directly overlay the transformed second-modality medical image with the first-modality medical image for intuitive observation of the local MRI location within the overall CT scan. Another format is to output the transformation parameters themselves, including the scale factor, rotation matrix, and translation vector, for use by other systems or subsequent processing. Yet another format is to output a transformed second-modality medical image mask or segmentation label, mapping it to the first-modality medical image space for fusion analysis or surgical planning.

[0097] Preferably, for the registration scenario of whole-body CT and local MRI of the lower extremities, the output result can be represented as the local MRI precisely covering the anatomical structures of the corresponding lower extremity. For example, when the second modal medical image is a left knee joint MRI, after registration, this MRI image will be accurately placed in the joint region between the distal femur and proximal tibia on the left side of the whole-body CT, achieving spatial consistency between bony and soft tissue structures. Through the above output, the precise localization and alignment of the region of interest in the whole-body anatomical reference space is achieved, providing a foundation for subsequent multimodal image fusion analysis, lesion localization, and surgical navigation.

[0098] In a specific embodiment, the following discloses a specific embodiment of the manifold learning-based multimodal lower limb medical image registration method applicable to this application, which specifically includes the following steps: Let the first modality image covering the entire area of ​​both lower limbs (the CT image covering the entire area of ​​both lower limbs) be denoted as The second modality image (MRI image covering a local area) is denoted as Let voxels be used. The corresponding physical coordinates are ,in Indicates the left and right coordinates. Indicates the coordinates of the forward and backward directions. This represents the head-to-feet orientation coordinates. Let the voxel spacing between the two modes in the unified physical orientation coordinate system be... and Preferably, the two modes are first unified to the same physical direction coordinate system before subsequent processing.

[0099] Let the surface mask of the first modality image / overall image be denoted as . The target mask of the second modality image / local image is .in, and It can be obtained by body surface threshold segmentation and significant connected component screening.

[0100] For the left-right separation of the overall image, if There are two main connected components. and Then the physical centroid of each connected domain is defined as equation (1): ; Compare the left and right coordinates of the centroids of two connected regions, and select the one with the larger left-hand coordinate as the candidate region. The region with the smaller horizontal coordinate is the candidate region on the right. .

[0101] like If it only represents a single primary connected component, then first construct its two-dimensional projection point set. Then divide it into two cluster subsets. and Based on this, define the background dividing lines. The optimization objective is given by equation (2): ; in, This represents the cost function for the dividing line passing through the background region. Represents the smoothing constraint coefficient. This represents the value of the dividing polyline at position y. By minimizing equation (2), we obtain a dividing polyline that passes through the gap in the real background between the two lower limbs as much as possible. Figure 2 As an example, the middle section shows a dividing line that separates the lower limbs into two cluster subsets. and .

[0102] To reduce the impact of high spatial resolution heterogeneity on feature comparability, the voxel spacing of the overall image is used. As a reference, a local image is resampled to obtain a uniform local image. It can be expressed as equation (3): ; in, This represents a 3D resampling operator based on the target voxel spacing. Preferably, it assumes a local image. The original voxel spacing is The original voxel size is The target voxel size after resampling Based on the principle of maintaining consistency in physical scope, it is determined as follows: ; Subsequently, under a unified physical coordinate system, A resampling grid is constructed using the target voxel spacing, and three-dimensional linear interpolation is used to calculate the gray value at each target voxel, thereby obtaining a uniform local image. For target voxels that are outside the domain of the original local image, their grayscale values ​​are preferably assigned to a preset background value.

[0103] After the candidate regions are determined, high-dimensional features of local blocks are extracted from the corresponding candidate regions of the overall image and the target regions of the local image. Let... Indicated by voxels Centered on, side length determined by parameters The controlled three-dimensional neighborhood block, then the local block feature vector It can be expressed as equation (5): in, This represents expanding a 3D neighborhood block into a high-dimensional vector. The high-dimensional local block set is composed of the feature vectors corresponding to all valid voxels. .

[0104] Construct a symmetric nearest neighbor graph on a high-dimensional local block set. For any two sampling points... and The graph edge weights are defined by equation (6): ; when belong The condition holds true only if the nearest neighbor set is used; otherwise... .in, and These represent the sampling points. and sampling points The local scale parameter is defined as follows: / The average distance between it and all its nearest neighbors. This parameter allows the similarity measure to adapt to changes in the density of the local point set.

[0105] Further define the weighted adjacency matrix. The corresponding degree matrix The diagonal elements satisfy Off-diagonal elements satisfy Based on this, construct the Graph Laplace matrix. Its satisfaction By solving the generalized characteristic equation (7): ; The low-dimensional embedding coordinates of the high-dimensional local block point set can be obtained. Among them, Let m be the generalized eigenvalue corresponding to the eigenvector m. Preferably, eigenvectors corresponding to several non-trivial eigenvalues ​​are selected. As a manifold embedding basis, the feature vector components are mapped back to their corresponding voxel positions according to the positions of each sampling point in the original voxel space, thereby forming a multi-channel manifold feature map. and ,in Manifold feature map representing the overall candidate region. A manifold feature map representing a local target region.

[0106] Since the manifold embeddings of different modal images are obtained independently, their low-dimensional feature representations are not unique. That is, after independently solving the manifold embeddings, the low-dimensional feature representations obtained for different modal images lack cross-modal consistency. Specifically, on the one hand, the feature vector corresponding to any embedding channel can be multiplied as a whole without changing the manifold structure. This leads to ambiguity in channel notation; on the other hand, there may be differences in orthogonal bases such as rotation or reflection between low-dimensional embedding bases obtained from different modes, causing inconsistent directional representations of the same structural mode in the low-dimensional coordinate system. To eliminate these differences, let the embedding matrices corresponding to the global candidate region and the local target region be respectively... and First, determine the symbol correction vector based on channel correlation. ,in Then solve for the orthogonal alignment matrix. As shown in equation (8): ; in, To meet The orthogonality matrix of the constraints is solved during the optimization process. The optimal orthogonal alignment matrix is ​​used to obtain the manifold feature representation of the aligned local image. , Represented by vector The elements are a diagonal matrix composed of diagonal elements. It should be noted that this step completes the representation unification in the feature space of the multimodal manifold, rather than the registration result in the final physical space.

[0107] like Figures 3a to 3f As shown, Figures 3a to 3c (CT candidate feature channels 1–3) correspond to the intensity distribution of the first three nontrivial feature channels obtained by Graph Laplacian embedding in the candidate region of the left or right lower limb of the first modality image (CT); Figures 3e to 3f (MRI Alignment Feature Channels 1–3) correspond to the second modality image (MRI) generated in the target anatomical region (ipsilateral lower limb) through the same manifold learning process.

[0108] After manifold alignment, let the overall candidate region manifold feature map be... The manifold feature map after local target region alignment is as follows The coarse transform from local image to global candidate region is estimated using Fourier-Merlin transform. Its formal definition is Equation (9): in, Indicates the scale factor. Represents the rotation matrix. This represents the translation vector. Obtained by matching the radial distribution of the spectrum. Obtained through multi-level angle search. The candidate regions are determined jointly by phase correlation and translation. This completes the coarse localization of the local image relative to the global candidate region.

[0109] To further improve registration accuracy based on coarse localization, a distance map can be constructed for the corresponding candidate regions. and in coarse transformation Nearby restricted neighborhood Solve the local refinement transformation. Define the local refinement objective function. For equation (10): ; in, The target mask region of the local image to be registered. To apply the transformation, Indicates the physical centroid of the corresponding region. The central constraint weights are represented. The final refinement transformation is determined by equation (11): ; Among them, the restricted neighborhood Constraints are imposed by translation amplitude thresholds, local attitude perturbation thresholds, and allowable center drift thresholds to avoid re-entering a large-scale global search. Only when... Only when the refinement result remains within a reasonable overlap range of the current candidate lateral regions will the refinement result be accepted as the final spatial transformation.

[0110] Ultimately, it will transform. Applying this to the original local image, the medical image registration result in the global image coordinate system is obtained, as shown in Equation (12): ; in, This represents a local image after medical image registration. For lower limb scenes, the registered local image or its corresponding label, mask, and other derived results can be backfilled into the overall image coordinate space to obtain the final multimodal fusion output. Figures 4a to 4c (Slice views in axial, coronal, and sagittal planes, respectively) are an example of the registration results. The bright white / grayish dense areas are soft tissues as shown on MRI, while the medium-low gray non-dense areas are bony structures of the entire lower limb as shown on CT.

[0111] Through the above-mentioned variable definition, lateral determination, manifold construction, feature alignment, coarse localization and local refinement process, the present invention can achieve multimodal medical image registration for lower limb scenarios under conditions of inconsistent anatomical coverage, highly heterogeneous spatial resolution and large differences in body position.

[0112] Further reference Figure 5 As an implementation of the above-described method, this application provides an embodiment of a multimodal lower limb medical image registration system based on manifold learning. This system embodiment is similar to... Figure 1Corresponding to the method embodiments shown, the system can be specifically applied to various electronic devices.

[0113] refer to Figure 5 A multimodal lower limb medical image registration system based on manifold learning, comprising: The medical image acquisition module 201 is configured to acquire a first modal medical image and a second modal medical image to be registered. The first modal medical image is a CT image covering the entire range of both lower limbs, and the second modal medical image is an MRI image covering a local area. The lateral identification module 202 is configured to extract local limb masks from the second modality medical image and determine its lateral identification, perform left-right separation on the first modality medical image, and obtain the corresponding side candidate region corresponding to the lateral identification of the second modality medical image; The resampling module 203 is configured to resample the second modality medical image with reference to the spatial resolution of the first modality medical image, and perform spatial resolution uniformity and preprocessing. The manifold feature map acquisition module 204 is configured to extract high-dimensional feature vectors of local three-dimensional neighborhood blocks centered on voxels in the corresponding candidate region and the region corresponding to the local limb mask, respectively, and construct a high-dimensional local block feature point set; construct a symmetric nearest neighbor graph based on the high-dimensional local block feature point set and determine the edge weights to form a weighted adjacency graph; construct a graph Laplacian matrix based on the weighted adjacency graph, solve for the low-dimensional embedding coordinates, and map them back to the original spatial position to obtain a first multi-channel manifold feature map and a second multi-channel manifold feature map; The manifold alignment module 205 is configured to perform manifold alignment on the first multi-channel manifold feature map and the second multi-channel manifold feature map to eliminate the representation ambiguity caused by independent embedding of different modalities and obtain a set of aligned manifold feature maps. The coarse registration transformation module 206 is configured to input the aligned set of manifold feature maps into the frequency domain medical image registration module, and estimate the scale factor, rotation parameter and translation parameter based on the Fourier-Merlin transform to obtain the coarse medical image registration transformation; The fine registration transformation module 207 is configured to perform restricted neighborhood search and refinement optimization on the overlapping area of ​​the region corresponding to the local limb mask and the corresponding side candidate region, based on the coarse medical image registration transformation, to obtain the final spatial transformation. The registration module 208 is configured to apply the final spatial transformation to the second modality medical image, map it to the spatial coordinate system of the first modality medical image, and output the registration result.

[0114] In another aspect, this application also provides a computer-readable storage medium, which may be included in the electronic device described in the above embodiments; or it may exist independently and not assembled into the electronic device. The aforementioned computer-readable storage medium carries one or more programs, which, when executed by the electronic device, cause the electronic device to perform the following... Figure 1 The method shown.

[0115] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0116] Through the above description of the embodiments, those skilled in the art can clearly understand that the methods of the above embodiments can be implemented by means of software plus necessary general-purpose hardware platforms. Of course, they can also be implemented by hardware, but in many cases the former is a better implementation method. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product is stored in a storage medium (such as ROM / RAM, magnetic disk, optical disk) and includes several instructions to cause a computer terminal (which may be a mobile phone, computer, server, or network device, etc.) to execute the methods described in the various embodiments of the present invention.

[0117] The embodiments of the present invention have been described above with reference to the accompanying drawings. However, the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of the present invention without departing from the spirit and scope of the claims. All of these forms are within the protection scope of the present invention.

Claims

1. A multimodal lower limb medical image registration method based on manifold learning, characterized in that, Includes the following steps: S1: Acquire the first and second modal medical images to be registered. The first modal medical image is a CT image covering the entire range of both lower limbs, and the second modal medical image is an MRI image covering a local area. S2: Extract local limb masks from the second modality medical image and determine its side orientation; perform left-right separation on the first modality medical image and obtain the corresponding side candidate region corresponding to the side orientation of the second modality medical image; S3: Using the spatial resolution of the first modality medical image as a reference, resample the second modality medical image, and perform spatial resolution uniformity and preprocessing; S4: Extract high-dimensional feature vectors of local 3D neighborhood blocks centered on voxels from the corresponding candidate regions and the regions corresponding to the local limb masks, respectively, and construct a high-dimensional local block feature point set; construct a symmetric nearest neighbor graph based on the high-dimensional local block feature point set and determine the edge weights to form a weighted adjacency graph; construct a graph Laplacian matrix based on the weighted adjacency graph, solve for the low-dimensional embedding coordinates, and map them back to the original spatial position to obtain the first multi-channel manifold feature map and the second multi-channel manifold feature map; S5: Perform manifold alignment on the first multi-channel manifold feature map and the second multi-channel manifold feature map to eliminate the representation ambiguity caused by independent embedding of different modalities and obtain a set of aligned manifold feature maps; S6: Input the aligned set of manifold feature maps into the frequency domain medical image registration module, and estimate the scale factor, rotation parameter and translation parameter based on Fourier-Merlin transform to obtain the coarse medical image registration transformation; S7: Based on the coarse medical image registration transformation, a restricted neighborhood search and refinement optimization are performed on the overlapping area between the region corresponding to the local limb mask and the corresponding side candidate region to obtain the final spatial transformation. S8: Apply the final spatial transformation to the second modality medical image, map it to the spatial coordinate system of the first modality medical image, and output the registration result.

2. The multimodal lower limb medical image registration method based on manifold learning according to claim 1, characterized in that, In step S2, performing left-right separation on the first modality medical image includes: Extract the lower limb surface mask from the first modality of medical images; If the lower limb surface mask contains two main connected regions, then the left and right are marked according to the horizontal coordinate of the physical centroid. The candidate region with the larger horizontal coordinate is the left candidate region, and the candidate region with the smaller horizontal coordinate is the right candidate region. If the lower limb surface mask is a single connected region, then perform bi-cluster clustering on the two-dimensional projection point set of the lower limb surface mask, estimate the separating polylines in the background gaps between the bi-clusters, and perform separation based on the separating polylines.

3. The multimodal lower limb medical image registration method based on manifold learning according to claim 1, characterized in that, In step S2, extracting a local limb mask from the second modality medical image and determining its side includes: Extract local limb masks from the second modality of medical images; Based on the horizontal coordinates of the physical centroid of the connected domain of the local limb mask or the morphological matching score of the local limb mask, the second modality medical image is determined to be left, right, or bilateral. The morphological matching score is calculated as a weighted sum of the following indicators based on the local limb mask: The similarity of shapes, consistency of physical dimensions, consistency of volume ratios, and consistency of center position after alignment in the main direction.

4. The multimodal lower limb medical image registration method based on manifold learning according to claim 1, characterized in that, In step S3, the resampling uses three-dimensional linear interpolation to construct a resampling grid based on the voxel spacing of the first modality medical image, and constrains the size of the target voxel after resampling based on the principle of maintaining physical range consistency.

5. The multimodal lower limb medical image registration method based on manifold learning according to claim 1, characterized in that, In step S4, the edge weights of the nearest neighbor graph are defined as follows: When two sampling points are neighbors, the edge weight is equal to the value of an exponential function with the square of the difference between the eigenvectors of the two points divided by the negative of the product of the local scale parameters of the two points; otherwise, the edge weight is zero. The local scale parameter of each sampling point is defined as the average distance between the current point and all its neighbors.

6. The multimodal lower limb medical image registration method based on manifold learning according to claim 5, characterized in that, In step S5, performing manifold alignment on the first multi-channel manifold feature map and the second multi-channel manifold feature map includes: Let the embedding matrix of the corresponding candidate region be the first embedding matrix, and the embedding matrix of the region corresponding to the local limb mask be the second embedding matrix. First, determine the symbol correction vector based on the correlation between each channel. Each component of the symbol correction vector is positive or negative. Then, solve for an orthogonal alignment matrix such that the Frobenius norm of the difference between the first embedding matrix and the second embedding matrix multiplied by the symbol correction diagonal matrix and then multiplied by the orthogonal alignment matrix is ​​minimized, thus obtaining a set of aligned manifold feature maps.

7. The multimodal lower limb medical image registration method based on manifold learning according to claim 6, characterized in that, The aligned set of manifold feature maps includes a first multi-channel manifold feature map that retains the original manifold embedding, and an aligned manifold feature map after performing a manifold alignment transformation on the second multi-channel manifold feature map; The aligned manifold feature map is obtained by multiplying the second embedding matrix by the sign-corrected diagonal matrix and the orthogonal alignment matrix in sequence.

8. The multimodal lower limb medical image registration method based on manifold learning according to claim 1, characterized in that, In step S6, based on the Fourier-Merlin transform, the scale factor, rotation parameters, and translation parameters are estimated to obtain the coarse medical image registration transform, including: The aligned set of manifold feature maps are unified into a common computational grid. The scaling factor is estimated based on the radial distribution of the spectrum, the rotation parameters are estimated based on the multi-level angle search of the spectrum amplitude, and the translation parameters are determined by combining phase correlation and translation candidate evaluation.

9. A multimodal lower limb medical image registration method based on manifold learning according to any one of claims 1-8, characterized in that, S7 includes: Obtain the body surface mask corresponding to the candidate region on the corresponding side and construct a distance map. Optimize within the preset translation and local pose perturbation range near the coarse medical image registration transformation. Only when the optimization result increases the overlap between the local limb mask and the body surface mask corresponding to the candidate region on the corresponding side or reduces the center offset, and the optimization displacement does not exceed the local drift threshold, is the optimization result accepted as the final spatial transformation.

10. A multimodal lower limb medical image registration system based on manifold learning, characterized in that, The system includes: The medical image acquisition module is configured to acquire a first modality medical image and a second modality medical image to be registered. The first modality medical image is a CT image covering the entire range of both lower limbs, and the second modality medical image is an MRI image covering a local area. The lateral identification module is configured to extract local limb masks from the second modality medical image and determine their laterality, perform left-right separation on the first modality medical image, and obtain the corresponding side candidate region corresponding to the laterality of the second modality medical image; The resampling module is configured to resample the second modality medical image with reference to the spatial resolution of the first modality medical image, and to perform spatial resolution uniformity and preprocessing. The manifold feature map acquisition module is configured to extract high-dimensional feature vectors of local three-dimensional neighborhood blocks centered on voxels from the corresponding candidate regions on the sides and the regions corresponding to the local limb masks, respectively, and construct a high-dimensional local block feature point set; construct a symmetric nearest neighbor graph based on the high-dimensional local block feature point set and determine the edge weights to form a weighted adjacency graph; construct a graph Laplacian matrix based on the weighted adjacency graph, solve for the low-dimensional embedding coordinates, and map them back to the original spatial position to obtain a first multi-channel manifold feature map and a second multi-channel manifold feature map; The manifold alignment module is configured to perform manifold alignment on the first multi-channel manifold feature map and the second multi-channel manifold feature map to eliminate representational ambiguity caused by independent embedding of different modalities and obtain a set of aligned manifold feature maps. The coarse registration transformation module is configured to input the aligned set of manifold feature maps into the frequency domain medical image registration module, and estimate the scale factor, rotation parameters and translation parameters based on the Fourier-Merlin transform to obtain the coarse medical image registration transformation; The fine registration transformation module is configured to perform restricted neighborhood search and refinement optimization on the overlapping area of ​​the region corresponding to the local limb mask and the corresponding side candidate region, based on the coarse medical image registration transformation, to obtain the final spatial transformation. The registration module is configured to apply the final spatial transformation to the second modality medical image, map it to the spatial coordinate system of the first modality medical image, and output the registration result.

Images