Non-rigid registration method of intraoperative ultrasound and preoperative image based on topological high-order mrf

Abstract

The application discloses a kind of based on topological high-order MRF intraoperative ultrasound and preoperative image non-rigid registration method, first iUS as target image T, preoperative image is as floating image M and is respectively handled;Then according to the Markov property of MRF, first define neighborhood system N and subgroup C, then define the context constraint of the target image T and floating image M after processing respectively by defined neighborhood system N and subgroup C;Again based on defined neighborhood system N and subgroup C, realize the non-rigid registration energy function of target image T and floating image M based on topological high-order MRF;Finally, based on the calculation of the non-rigid registration energy function of target image T and floating image M based on topological high-order MRF, based on multi-resolution additional variable reduction energy optimization algorithm is realized.The application improves the intraoperative dynamic decision-making ability of surgical navigation, provides accurate stereoscopic visual surgery guidance for physicians in clinic, and improves surgical precision.

Description

Technical Field

[0001] This invention relates to the fields of medical image navigation and image registration technology in medical image processing, specifically a non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF. Background Technology

[0002] Computer-assisted surgery utilizes computer technology to assist physicians in preoperative surgical planning, intraoperative surgical navigation, and related adjuvant therapies. It has been widely applied in clinical treatments such as neurosurgery, orthopedic surgery, spinal surgery, plastic surgery, and interventional oncology. Precise surgical navigation systems can significantly shorten surgical time and greatly improve surgical safety, serving as a crucial tool for achieving precise clinical diagnosis and treatment. They provide physicians with convenient, intelligent, and reliable surgical guidance in clinical practice.

[0003] Most current surgical navigation systems rely on preoperative images and lack the ability to make dynamic decisions during surgery. They cannot reflect changes in the structure and location of diseased tissues in a timely manner, leading to the failure of preoperative surgical planning and loss of basis for surgical implementation. Flexible organs such as brain tumors, liver, and kidneys are deformed and displaced by changes in intracranial pressure, respiration, heartbeat, collisions with surgical instruments, and traction, resulting in differences between the shape and position of organs during surgery and those in preoperative images. Although rigid transformations between organs during surgery and before surgery are compensated for through spatial orientation calibration, how to quantify the non-rigid transformations caused by deformation and displacement to provide surgeons with precise three-dimensional visualization for surgical guidance is a key technology that needs to be overcome in the future.

[0004] Thanks to the advantages of portable, inexpensive, and real-time ultrasound imaging, intraoperative ultrasound (iUS) has been widely used in recent years to acquire information on soft tissue deformation and drift during surgery. However, limited by the low contrast and high noise image quality of ultrasound images, iUS cannot be directly used for intraoperative surgical navigation. Therefore, a method to compensate for intraoperative soft tissue deformation and drift information by non-rigid registration of iUS with preoperative images has great research potential. With the continuous development of computer vision technology and theory, more accurate and efficient registration methods are becoming possible.

[0005] Image registration based on Markov Random Field (MRF) models is a non-rigid registration method using discrete optimization, requiring the selection of appropriate discrete optimization algorithms based on different forms of energy functions. The deformation model of MRF-based image registration algorithms is relatively simple, with FFD being the primary focus. The objective function and optimization algorithms are the main development directions for this type of algorithm. Current research mostly uses low-order MRFs with simple modeling, establishing correspondences through single-attribute matching in salient regions. While this achieves high computational efficiency, it suffers from significant errors when landmarks are unevenly distributed in the image space or when the registered images have large differences. High-order MRFs can more accurately describe complex systems and data structures, but their computational complexity is usually high, making precise reasoning difficult. Summary of the Invention

[0006] To address the shortcomings of existing technologies, the technical problem this invention aims to solve is to provide a non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF.

[0007] The technical solution of this invention to solve the aforementioned technical problem is to provide a non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF, characterized in that the method includes the following steps:

[0008] Step 1: Use iUS as the target image T and the preoperative image as the floating image M; then unify the size and resolution of the target image T and the floating image M, and then establish the target image T and the floating image M in the coordinate system of the target image T to obtain the processed target image T and the floating image M.

[0009] Step 2: Based on the Markov property of MRF, first define the neighborhood system N and sub-clique C, and then define the context constraints of the processed target image T and floating image M obtained in Step 1 through the defined neighborhood system N and sub-clique C respectively.

[0010] Step 3: Based on the neighborhood system N and sub-cluster C defined in Step 2, implement the non-rigid registration energy function for establishing the target image T and the floating image M based on the topological high-order MRF;

[0011] Step 4: Calculate the energy function of the non-rigid registration of the target image T and the floating image M based on the multi-resolution, variable-free reduced-order energy optimization algorithm.

[0012] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0013] (1) This invention can compensate for the soft tissue deformation and drift information during surgery to the preoperative image through non-rigid registration of iUS and preoperative image, thereby improving the dynamic decision-making ability of surgical navigation during surgery, solving the problem of image navigation failure based on preoperative image caused by tissue displacement during surgery, greatly reducing memory consumption, and providing doctors with precise three-dimensional visualization surgical guidance in clinical practice, thereby improving surgical accuracy.

[0014] (2) Based on the higher-order MRF theory, this invention establishes a non-rigid registration energy function between iUS and preoperative images through potential functions in neighborhood systems, subclusters and random fields.

[0015] (3) The unary clique potential function uses a set of multi-scale, multi-directional Gabor feature vectors to describe each voxel, capturing image texture information to improve registration accuracy. ZNCC is used as the similarity index for data items, giving full play to its anti-interference ability, and it can remain stable even in the presence of US image noise, thereby reducing registration error.

[0016] (4) By introducing the topological term of the higher-order group into the energy function, the complex relationship between the labels is described by the quaternary group, which constrains the smoothness and maintains the topological structure of the deformation field, so that the deformation is better constrained, the deformation field is accurately estimated, and the registration accuracy is improved.

[0017] (5) To address the difficulty in optimizing high-order MRF energy functions, this invention proposes a multi-resolution, variable-free energy optimization algorithm. This algorithm utilizes high-order variable terms that cannot become the global minimum and their corresponding values ​​as increments to the MRF potential function, thereby reducing the potential function's order without altering the variable values ​​corresponding to the potential function's minimum value. This method maintains topological invariance while ensuring image registration accuracy, thus achieving the solution for high-order MRF energy functions. Attached Figure Description

[0018] Figure 1 This is an overall flowchart of the present invention;

[0019] Figure 2 This is a diagram of the non-rigid registration algorithm of the present invention;

[0020] Figure 3 This is a neighborhood system diagram of the target image and floating image nodes in this invention;

[0021] Figure 4 This describes the structure of each sub-cluster in the two-dimensional and three-dimensional image neighborhood systems of the present invention;

[0022] Figure 5 This is a result of the image registration method in Embodiment 1 of the present invention. Detailed Implementation

[0023] Specific embodiments of the present invention are given below. These specific embodiments are only used to further illustrate the present invention in detail and do not limit the scope of protection of the present invention.

[0024] This invention provides a non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF (hereinafter referred to as the method), characterized by the following steps:

[0025] Step 1: Use iUS as the target image T and the preoperative image as the floating image M; then unify the size and resolution of the target image T and the floating image M, and then establish the target image T and the floating image M in the coordinate system of the target image T to obtain the processed target image T and the floating image M.

[0026] Preferably, in step 1, the preoperative image is a preoperative MRI image or a preoperative CT image.

[0027] Preferably, step 1 is performed in C3D software: the -resample-mm command is used to unify the size and resolution; the -reslice-identity command is used to reposition the floating image M space to the target image T space, initially establishing the transformation relationship between the two coordinate systems, providing good initial values ​​for non-rigid registration, and preventing it from getting trapped in local optima.

[0028] Step 2: Based on the Markov property of MRF, first define the neighborhood system N and sub-clique C, and then define the context constraints of the processed target image T and floating image M obtained in Step 1 through the defined neighborhood system N and sub-clique C respectively.

[0029] Preferably, in step 2, in image processing, the Markov property refers to the fact that pixels or feature points in an image are only affected by their spatially neighboring points and are independent of other points in the image.

[0030] Preferably, in step 2, contextual constraints are a method to enhance image analysis and processing capabilities by utilizing the relationships between pixels or regions in an image. Pixels or regions in an image do not exist in isolation but are interconnected. This interconnectivity can be utilized through contextual constraints to help better understand and process the image.

[0031] Preferably, in step 2, the neighborhood system N defines the neighborhood relationships of each pixel in the image to determine which pixels have direct dependencies, thereby introducing local context information into the model; the subclique C is a set of pixels in the image used to define potential functions that describe the interactions between pixels. In MRF, potential functions are typically defined on the subclique C to capture the dependencies between pixels. Based on the size of the clique (i.e., the number of points p contained in the clique), the unary clique C1, binary clique C2, and quaternary clique C4 are defined as follows:

[0032] C1={γ|γ∈Γ}

[0033] C2={(γ,η)|γ∈Γ,η∈N γ}

[0034] C4={(γ,η,μ,ψ)|γ∈Γ,η∈N γ ,μ∈N η ,ψ∈N μ ,η∈N μ ,μ∈N ψ} (1)

[0036] In equation (1), i (i = γ, η, μ, ψ) are the control points in the image, also known as voxels; N i The neighborhood system representing voxel i (i = γ, η, μ, ψ).

[0037] Step 3: Based on the neighborhood system N and sub-cluster C defined in Step 2, implement the non-rigid registration energy function E for establishing the target image T and the floating image M based on the topological high-order MRF, i.e., the MRF energy function E. MRF ;

[0038] Preferably, in step 3, the minimum value of the non-rigid registration energy function E(T) is sought using a non-rigid registration algorithm, that is, to find the minimum value of the voxels n∈Ω in the target image. t Mapping to floating image T(n)∈Ω m The transformation Trans is an approximation of the non-rigid registration energy function in the MRF, i.e., the MRF energy function E. MRF :

[0039]

[0040] In equation (2), the MRF energy function E MRF P represents the data terms (i.e., the potential function of the unary group C1) for each voxel i (i = γ, η, μ, ψ). γ The smoothing term (i.e., the potential function of the binary clique C2) P γη And the topological term (i.e., the potential function of higher-order cliques) P γημψ The system is composed of two parts; where the set of deformation displacement vectors for each voxel in the target image T and the floating image M is defined as L = {l γ ,l η ,l μ ,…,l ψ The potential function (}) is the unit of search and optimization in the MRF method. It can capture complex relationships between multiple pixels, further enhancing contextual constraints.

[0041] Preferably, in step 3, a set of multi-scale, multi-directional Gabor feature vectors A(·) are used to describe each voxel in the data item to capture image texture information and improve registration accuracy. For multi-scale, high-frequency features are set to obtain global fine-scale edges, and low-frequency features are set to obtain local coarse-scale edges. The parameters of the high-frequency features are: wavelength 5-10, frequency 0.1-0.2, and bandwidth 0.05, and the parameters of the low-frequency features are: wavelength 10-20, frequency 0.05-0.1, and bandwidth 0.4. For multi-directional, the direction range of the Gabor filter is transformed from 0, π / 4, π / 2, 3π / 4 to π to obtain image edge information from vertical to horizontal.

[0042]

[0043] In equation (3), η -1 (·) represents the mapping of node γ∈Γ to each ordinary point n∈Ω. t The inverse mapping function; ZNCC is used as the similarity index for data items because it has strong anti-interference properties and remains stable even in the presence of noise in the US image, thereby reducing registration errors. Figure 3 As shown, the registration performance is measured by the ratio of the adjacent neighborhood (AN) and the boundary neighborhood (BN) of node u in the target image T and node T(u) in the floating image M. The higher the average similarity of the adjacent neighborhood and the lower the average similarity of the boundary neighborhood, the more reliable the matching is.

[0044] Preferably, in step 3, in the smoothing term, to ensure the deformation field satisfies a certain smoothness, an intuitive prior condition is adopted: the labels corresponding to two adjacent nodes in the neighborhood should be as consistent as possible; based on this prior condition, the binary clique C2 is used to define the potential function that applies smoothing constraints to the deformation field:

[0045] P γη (l γ ,l η )=min(λ‖l γ -lx‖2,c) (4)

[0046] The minimum value of the binary clique potential function is defined as the L2 norm of the vector difference of the binary clique, in equation (4). γ ,l η Let γ and η represent the labels of two adjacent nodes γ and η in the target image T, respectively. λ = 2 represents the smoothness weight, and c = 0.1 is used to control the maximum cost.

[0047] Preferably, in step 3, in the topology term, the values ​​of the angular Jacobian determinants corresponding to the four vertices of the basic unit, i.e., the quaternary clique C4, are all greater than 0 to maintain the topological structure of the deformation field; after expanding from two dimensions to three dimensions, the spatial topological structure constraint is expanded from the quaternary clique to the octet constraint, and the topological constraint condition is expanded from 4 dimensions to 64 dimensions; for thousands of nodes in the discrete domain, practical calculation is impractical. Preferably, as follows... Figure 4 As shown, to reduce the computational load, it is proposed to use 8 angular Jacobian determinants corresponding to 8 quaternions to approximate the complete 3D topology-preserving condition.

[0048] Preferably, in step 3, to avoid excessive cost from topological terms composed of higher-order cliques, a logarithmic function is proposed to penalize large deformations; for 3D images, the topological terms can be written as:

[0049]

[0050] In equation (5), w is the cost value when the topology is not retained, satisfying log(J(l γ ,l η ,l μ ,l ψ )+1)<<w<1.

[0051] Step 4: Calculate the energy function of the non-rigid registration of the target image T and the floating image M based on the multi-resolution, variable-free reduced-order energy optimization algorithm.

[0052] Preferably, step 4 specifically involves: using higher-order variable terms that cannot become the global minimum and their corresponding values ​​as increments to the MRF potential function; these increments reduce the order of the potential function without changing the variable values ​​corresponding to the minimum value of the potential function; these increments are expressed as polynomials:

[0053]

[0054] In equation (6), C h k represents a higher-order clique. C For real coefficients, we use an exhaustive search method to find higher-order variable terms ω for the MRF potential function.

[0055] Preferably, in step 4, two Intel Xeon Gold 6226R (16 cores / 16 threads) CPUs with a dual-socket configuration are used for optimized calculations.

[0056] Preferably, in step 4, a multi-resolution deformation field sampling strategy is adopted to achieve non-rigid MRF registration; the input image is processed at ratios of 1 / 4, 1 / 2, and 1 / 1, using three levels of processing to achieve sparse sampling from coarse label distribution in a large search space to dense sampling from fine label distribution in a small search space; low-level registration obtains the deformation field based on the energy function; intermediate and high-level registration obtain the deformation field through the energy function and previous deformations, and upsample the multisaliency map calculated at the previous resolution; the multi-resolution processing sampling strategy transitions from sparse sampling in a wide search region to dense sampling in a limited region, thereby optimizing registration accuracy and computational efficiency.

[0057] Example 1:

[0058] For cases 5, 12, and 17 in the Resect public dataset, the specific steps are as follows:

[0059] In step 1, the target image T and the floating image M are 256×256×256 in size and 0.5mm×0.5mm×0.5mm in resolution.

[0060] In step 3, the parameters for the high-frequency features are wavelength 5, frequency 0.2, and bandwidth 0.05. The parameters for the low-frequency features are wavelength 10, frequency 0.1, and bandwidth 0.4.

[0061] Figure 5 In the dataset, a1 represents the pre-registration MRI and iUS images of case 5 in the RESECT dataset; a2 represents the post-registration MRI and iUS images of case 5 in the RESECT dataset; b1 represents the pre-registration MRI and iUS images of case 12 in the RESECT dataset; b2 represents the post-registration MRI and iUS images of case 12 in the RESECT dataset; c1 represents the pre-registration MRI and iUS images of case 17 in the RESECT dataset; and c2 represents the post-registration MRI and iUS images of case 17 in the RESECT dataset. Blue arrows indicate the boundaries of lesions in MRI images, and red arrows indicate the boundaries of lesions in iUS images.

[0062] Depend on Figure 5 As can be seen, the boundary of the lesion in the pre-registration MRI image was not aligned with the boundary of the lesion in the iUS image. After registration, the boundary of the lesion in the MRI image changed non-linearly, achieving alignment with the boundary of the lesion in the iUS image. Any aspects not described in this invention are applicable to existing technologies.

Claims

1. A non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF, characterized in that, The method includes the following steps: Step 1: Use iUS as the target image T and the preoperative image as the floating image M; then unify the size and resolution of the target image T and the floating image M, and then establish the target image T and the floating image M in the coordinate system of the target image T to obtain the processed target image T and the floating image M. Step 2: Based on the Markov property of MRF, first define the neighborhood system N and sub-clique C, and then define the context constraints of the processed target image T and floating image M obtained in Step 1 through the defined neighborhood system N and sub-clique C respectively. Step 3: Based on the neighborhood system N and sub-cluster C defined in Step 2, implement the non-rigid registration energy function for establishing the target image T and the floating image M based on the topological high-order MRF; Seeking to extract voxels from the target image Mapping to floating image The transformation Trans is an approximation of the non-rigid registration energy function in the MRF, i.e., the MRF energy function. : （2） In equation (2), the MRF energy function It is each voxel Data items Smooth items and topology terms The system is composed of two parts; where the set of deformation displacement vectors for each voxel in the target image T and the floating image M is defined as follows: , is the unit of search and optimization in the MRF method; ; In the data items, a set of multi-scale, multi-directional Gabor feature vectors are used. Each voxel is described to capture image texture information to improve registration accuracy. For multiple scales, high-frequency features are set to obtain global fine-scale edges, and low-frequency features are set to obtain local coarse-scale edges. The parameters of the high-frequency features are: wavelength 5~10, frequency 0.1~0.2, bandwidth 0.05, and the parameters of the low-frequency features are: wavelength 10~20, frequency 0.05~0.1, bandwidth 0.

4. For multiple directions, the direction range of the Gabor filter is transformed from 0, π / 4, π / 2, 3π / 4 to π to obtain image edge information from vertical to horizontal. （3） In equation (3), ; Representative node Mapped to each ordinary point The inverse mapping function; AN represents the nearest neighbor; BN represents the boundary neighbor; Step 4: Calculate the energy function of the non-rigid registration of the target image T and the floating image M based on the multi-resolution, variable-free reduced-order energy optimization algorithm.

2. The non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF according to claim 1, characterized in that, In step 1, the preoperative images are either MRI or CT images.

3. The non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF according to claim 1, characterized in that, Step 1 is performed in C3D software: The -resample-mm command is used to unify the size and resolution; the -reslice-identity command is used to reposition the floating image in M ​​space to the target image in T space, initially establishing the transformation relationship between the two coordinate systems, providing good initial values ​​for non-rigid registration, and preventing it from getting trapped in local optima.

4. The non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF according to claim 1, characterized in that, In step 2, within the MRF, based on the size of the clique, unary clique C1, binary clique C2, and quaternary clique C4 are defined as follows: （1） In equation (1), Points that control an image are called voxels; Representative voxels The neighborhood system; .

5. The non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF according to claim 1, characterized in that, In step 3, an intuitive prior condition is adopted in the smoothing term: the labels corresponding to two adjacent nodes in the neighborhood should be as consistent as possible; based on this prior condition, the binary clique C2 is used to define the potential function that applies smoothing constraints to the deformation field. （4） The minimum value of the binary clique potential function is defined as the L2 norm of the vector difference of the binary clique, in equation (4). These represent two adjacent nodes in the target image T. and The tag, =2 represents the smoothness weight, and c=0.1 is used to control the maximum cost.

6. The non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF according to claim 1, characterized in that, In step 3, in the topology term, the values ​​of the angular Jacobian determinants corresponding to the four vertices of the quaternary clique C4 in the two-dimensional image are all greater than 0 to maintain the topology of the deformation field; after expanding from two dimensions to three dimensions, the spatial topology constraint is expanded from the quaternary clique to the octet constraint, and the topology constraint condition is expanded from 4 dimensions to 64 dimensions. In step 3, to avoid excessive cost from topological terms composed of higher-order cliques, a logarithmic function is proposed to penalize large deformations; for 3D images, the topological terms can be written as: （5） In equation (5), w is the cost value when the topology is not retained, which satisfies .

7. The non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF according to claim 1, characterized in that, Step 4 specifically involves using higher-order variable terms that cannot become the global minimum and their corresponding values ​​as increments to the MRF potential function. These increments reduce the order of the potential function without changing the variable values ​​corresponding to the minimum value of the potential function. This increment is expressed as a polynomial: （6） In equation (6), Indicates a higher-order group. Given real coefficients, an exhaustive search method is used to find higher-order variable terms for the MRF potential function. .

8. The non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF according to claim 1, characterized in that, In step 4, a multi-resolution deformation field sampling strategy was adopted to achieve non-rigid MRF registration; the input image was processed at ratios of 1 / 4, 1 / 2 and 1 / 1, and a three-level processing was adopted to achieve sparse sampling from coarse label distribution in a large search space to dense sampling from fine label distribution in a small search space. Low-level registration obtains the deformation field based on the energy function; intermediate and high-level registration obtain the deformation field through the energy function and previous deformation, and upsample the multisaliency map calculated from the previous resolution. The multi-resolution processing sampling strategy transitions from sparse sampling over a wide search area to dense sampling over a limited area, thereby optimizing registration accuracy and computational efficiency.

Images