Intraoperative ultrasound and preoperative image non-rigid registration method based on topological high-order MRF (Markov Random Field)

Abstract

The invention discloses an intraoperative ultrasound and preoperative image non-rigid registration method based on topological high-order MRF (Markov Random Field). The method comprises the following steps: firstly, taking iUS as a target image T, taking a preoperative image as a floating image M, and respectively processing the target image T and the floating image M; then according to the Markov property of the MRF, firstly defining a neighborhood system N and a sub-group C, and then respectively defining context constraints of the processed target image T and the floating image M through the defined neighborhood system N and the defined sub-group C; establishing a non-rigid registration energy function of the target image T and the floating image M based on the topological high-order MRF on the basis of the defined neighborhood system N and the sub-group C; and finally, realizing calculation of a non-rigid registration energy function of the target image T and the floating image M based on the topological high-order MRF based on a multi-resolution additional-variable-free reduced-order energy optimization algorithm. According to the method, the intraoperative dynamic decision-making capability of surgical navigation is improved, precise three-dimensional visual surgical guidance is provided for doctors clinically, and the surgical precision is improved.

Description

Technical Field

[0001] The present invention relates to the field of medical image navigation and the technical field of image registration in medical image processing, and specifically to a non-rigid registration method between intraoperative ultrasound and preoperative images based on topological high-order MRF. Background Art

[0002] Computer-assisted surgery (CAD) utilizes computer technology to assist physicians with preoperative surgical planning, intraoperative surgical navigation, and related adjunctive treatments. It has been widely used in clinical procedures such as neurosurgery, orthopedics, spinal surgery, plastic surgery, and interventional tumor treatment. Precise surgical navigation systems can significantly shorten surgical time and significantly improve surgical safety. They are a crucial tool for achieving precise clinical diagnosis and treatment, providing physicians with convenient, intelligent, and reliable surgical guidance.

[0003] Most current surgical navigation systems, based on preoperative imaging, lack dynamic intraoperative decision-making capabilities and are unable to promptly reflect changes in the structure and position of lesions during surgery, leading to the failure of preoperative surgical planning and the loss of the basis for surgical implementation. Flexible tissues and organs such as brain tumors, livers, and kidneys deform and displace due to changes in intracranial pressure, respiration, heartbeat, collisions with surgical instruments, and traction, resulting in differences in the shape and position of the organs during surgery and those in preoperative imaging. Although the rigid transformation between the organs during and before surgery is compensated after spatial position calibration, quantifying the non-rigid transformations caused by deformation and displacement to provide doctors with precise three-dimensional visual surgical guidance is a key technology that needs to be overcome in the future.

[0004] In recent years, intraoperative ultrasound (iUS) has been widely used to obtain intraoperative soft tissue deformation and drift information, thanks to its portability, low cost, and real-time imaging. However, due to the low contrast and high noise quality of ultrasound images, iUS cannot be directly used for intraoperative surgical navigation. Therefore, a method for compensating intraoperative soft tissue deformation and drift information to preoperative images through non-rigid registration of iUS with preoperative images has great research prospects. With the continuous development of computer vision technology and theory, more accurate and efficient registration methods are becoming possible.

[0005] Image registration based on the Markov Random Field (MRF) model is a discrete optimization, non-rigid registration method that requires selecting an appropriate discrete optimization algorithm based on different forms of energy functions. MRF-based image registration algorithms have a relatively simple deformation model, FFD, and the objective function and optimization algorithm are the main development directions of this type of algorithm. Most current studies use simple low-order MRFs to establish correspondences by matching single attributes in salient regions. While this can achieve high computational efficiency, it suffers from large errors when landmarks are unevenly distributed in the image space and when the registered images differ significantly. Higher-order MRFs can more accurately describe complex systems and data structures, but their computational complexity is generally high, making precise reasoning difficult. Summary of the Invention

[0006] In view of the shortcomings of the existing technology, the technical problem to be solved by the present invention is to provide a non-rigid registration method between intraoperative ultrasound and preoperative images based on topological high-order MRF.

[0007] The technical solution of the present invention to solve the 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 comprises the following steps:

[0008] Step 1: Use the 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 target image T coordinate system to obtain the processed target image T and the floating image M;

[0009] Step 2: According to the Markov property of MRF, first define the neighborhood system N and subcluster C, and then use the defined neighborhood system N and subcluster C to define the context constraints of the processed target image T and floating image M obtained in step 1 respectively;

[0010] Step 3: Based on the neighborhood system N and subcluster C defined in step 2, a non-rigid registration energy function is established between the target image T and the floating image M based on the topological high-order MRF.

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

[0012] Compared with the prior art, the present invention has the following beneficial effects:

[0013] (1) The present invention can compensate for the intraoperative soft tissue deformation drift information in the preoperative image through non-rigid registration of iUS and preoperative images, thereby improving the intraoperative dynamic decision-making ability of surgical navigation, solving the problem of image navigation failure based on preoperative images due to intraoperative tissue displacement, greatly reducing memory consumption, and providing physicians with accurate three-dimensional visual surgical guidance in clinical practice, thereby improving surgical accuracy.

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

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

[0016] (4) By introducing the topological term of high-order clusters into the energy function, a four-element cluster is used to describe the complex relationship between labels, so that the smoothness is constrained and the topological structure of the deformation field is maintained, so that the deformation is better constrained, the deformation field is accurately estimated, and the accuracy of the registration is improved.

[0017] (5) In view of the difficulty in optimizing high-order MRF energy functions, this paper proposes a multi-resolution, variable-free, reduced-order energy optimization algorithm. This algorithm uses high-order variables that cannot reach the global minimum and their corresponding values as the added value of the MRF potential function, thereby reducing the order of the potential function without changing the variable values corresponding to the minimum value of the potential function. This method can maintain the invariance of the topological structure while ensuring the accuracy of image registration, thus achieving the solution of the high-order MRF energy function. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 It is the overall flow chart of the present invention;

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

[0020] Figure 3 A neighborhood system diagram of the target image and floating image nodes of the present invention;

[0021] Figure 4 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 diagram of the image registration method according to Example 1 of the present invention. DETAILED DESCRIPTION

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

[0024] The present 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 in that the method comprises the following steps:

[0025] Step 1: Use the 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 target image T coordinate system to obtain the processed target image T and the floating image M;

[0026] Preferably, in step 1, the preoperative image is an MRI preoperative image or a CT preoperative 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, and a preliminary transformation relationship between the two coordinate systems is established to provide a good initial value for non-rigid registration and prevent it from falling into a local optimum.

[0028] Step 2: According to the Markov property of MRF, first define the neighborhood system N and subcluster C, and then use the defined neighborhood system N and subcluster C to define the context constraints of the processed target image T and floating image M obtained in step 1 respectively;

[0029] Preferably, in step 2, in image processing, the Markov property means that a pixel point or feature point in an image is only affected by its spatial neighborhood points and is independent of other points in the image.

[0030] Preferably, in step 2, contextual constraints are a method for enhancing 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 interrelated. This interrelationship can be exploited through contextual constraints to help better understand and process the image.

[0031] Preferably, in step 2, the domain system N defines the neighborhood relationship of each pixel in the image, which is used to determine which pixels have direct dependencies, thereby introducing local context information into the model; the subcluster C is a collection of pixels in the image, which is used to define potential functions that describe the interactions between pixels. In MRF, the potential function is usually defined on the subcluster C to capture the dependencies between pixels. According to the size of the cluster (i.e., the number of points p contained in the cluster), the unary cluster C1, the binary cluster C2 and the four-element cluster C4 are defined as:

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

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

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

[0036] In formula (1), i (i = γ, η, μ, ψ) is the control point in the image, also called voxel; N i Represents the neighborhood system of voxel i (i = γ, η, μ, ψ).

[0037] Step 3: Based on the neighborhood system N and subcluster C defined in step 2, a non-rigid registration energy function, namely, the MRF energy function E, is established between the target image T and the floating image M based on the topological high-order MRF. MRF ;

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

[0039]

[0040] In formula (2), the MRF energy function E MRF is the data item of each voxel i (i = γ, η, μ, ψ) (i.e., the potential function of the unary cluster C1) P γ , smooth term (i.e., potential function of binary group C2) P γη and the topological term (i.e., the potential function of the higher-order group) P γημψ The target image T and the floating image M in each voxel deformation displacement vector set is defined as L = {l γ ,l η ,l μ ,…,l ψ}, is the unit of search and optimization in the MRF method. The potential function can capture the complex relationship between multiple pixels and further enhance the context constraint.

[0041] Preferably, in step 3, in the data item, a set of multi-scale and multi-directional Gabor feature vectors A(·) are used to describe each voxel to capture image texture information to 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, wherein 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 multi-direction, the direction range of the Gabor filter is set to be transformed from 0, π / 4, π / 2, 3π / 4 to π, so as to obtain image edge information from vertical to horizontal:

[0042]

[0043] In formula (3), η -1 (·) represents the mapping of node γ∈Γ to each common point n∈Ω t The inverse mapping function of ZNCC is used as the similarity index of data items because it has strong anti-interference ability and can remain stable even in the presence of US image noise, thereby reducing the registration error. Figure 3 As shown in Figure 1, the registration performance is measured by the ratio of the adjacent neighborhood (AN) and the boundary neighborhood (BN) of the node u in the target image T and the 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 match.

[0044] Preferably, in step 3, in the smoothness term, in order to make the deformation field satisfy a certain smoothness, an intuitive prior condition is adopted: the labels corresponding to two adjacent nodes in the neighborhood are as consistent as possible; based on this prior condition, the binary clique C2 is used to define the potential function that imposes a smooth constraint on the deformation field:

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

[0046] The minimum value of the binary cluster potential function is defined as the bi-norm of the vector difference of the binary cluster. In formula (4), l γ ,l η They represent the labels of two adjacent nodes γ and η in the target image T, λ = 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 in the basic unit, namely the quaternion C4, are all greater than 0 in the two-dimensional image to maintain the topological structure of the deformation field; after expanding from two dimensions to three dimensions, the spatial topological structure constraints are expanded from the quaternion to the octet constraint, and the topological constraints are expanded from 4 dimensions to 64 dimensions; for thousands of nodes in the discrete domain, actual calculation is unrealistic. Preferably, as Figure 4 As shown in the figure, in order to reduce the amount of calculation, it is proposed to use 8 angular Jacobian determinants corresponding to 8 quaternion groups to approximate the complete three-dimensional topology preservation conditions.

[0048] Preferably, in step 3, in order to avoid the high cost of the topological term composed of high-order groups, a logarithmic function is proposed to penalize large deformations; for a three-dimensional image, the topological term can be written as:

[0049]

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

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

[0052] Preferably, step 4 specifically comprises: using high-order variable terms and their corresponding assignments that cannot become global minimum values as the added value of the MRF potential function; the added value can reduce the order of the potential function and does not change the variable value corresponding to the minimum value of the potential function; the added value is expressed as a polynomial:

[0053]

[0054] In formula (6), C h represents a high-order group, k C is a real coefficient, and an exhaustive search method is used to find the high-order variable term ω for the MRF potential function.

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

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

[0057] Example 1:

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

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

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

[0061] Figure 5 In the figure, a1 is the MRI preoperative image and iUS of case 5 in the RESECT dataset before registration, a2 is the MRI preoperative image and iUS of case 5 in the RESECT dataset after registration, b1 is the MRI preoperative image and iUS of case 12 in the RESECT dataset before registration, b2 is the MRI preoperative image and iUS of case 12 in the RESECT dataset after registration, c1 is the MRI preoperative image and iUS of case 17 in the RESECT dataset before registration, and c2 is the MRI preoperative image and iUS of case 17 in the RESECT dataset after registration. The blue arrow indicates the boundary of the lesion in the MRI image, and the red arrow indicates the boundary of the lesion in the iUS image.

[0062] Depend on Figure 5 It can be seen that the boundary of the lesion in the MRI image before registration is not aligned with the boundary of the lesion in the iUS image. After registration, the boundary of the lesion in the MRI image changes nonlinearly and is aligned with the boundary of the lesion in the iUS image. The points not described in the present invention are applicable to the prior art.

Claims

1. A non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF, characterized by: The method comprises the following steps: Step 1: Use the 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 target image T coordinate system to obtain the processed target image T and the floating image M; Step 2: According to the Markov property of MRF, first define the neighborhood system N and subcluster C, and then use the defined neighborhood system N and subcluster C to define the context constraints of the processed target image T and floating image M obtained in step 1 respectively; Step 3: Based on the neighborhood system N and subcluster C defined in step 2, a non-rigid registration energy function is established between the target image T and the floating image M based on the topological high-order MRF. Step 4: Based on the multi-resolution variable-free reduced-order energy optimization algorithm, the non-rigid registration energy function of the target image T and the floating image M based on the topological high-order MRF is calculated.

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

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 M space to the target image T space, and a preliminary transformation relationship between the two coordinate systems is established to provide a good initial value for non-rigid registration and prevent it from falling into a local optimum.

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, in MRF, according to the size of the group, the unary group C1, the binary group C2 and the four-member group C4 are defined as: C1={γ|γ∈Γ} C2={(γ,η)|γ∈Γ,η∈N γ } C4={(γ,η,μ,ψ)|γ∈Γ,η∈N γ ,μ∈N η ,ψ∈N μ ,η∈N μ ,μ∈N ψ } (1) In formula (1), i (i = γ, η, μ, ψ) is the controlled point in the image, called voxel; N i Represents the neighborhood system of voxel i (i = γ, η, μ, ψ).

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, we seek to convert the voxel n∈Ω in the target image t Mapping to floating image T(n)∈Ω m The transformation Trans is the approximate approximation of the non-rigid registration energy function in MRF, that is, the MRF energy function E MRF : In formula (2), the MRF energy function E MRF is the data item P of each voxel i (i = γ, η, μ, ψ) γ , smooth term P γη and the topological term P γημψ The target image T and the floating image M in each voxel deformation displacement vector set is defined as L = {l γ ,l η ,l μ ,…,l ψ }, is the unit of search and optimization in the MRF method.

6. The non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF according to claim 5, characterized in that: In step 3, in the data item, a set of multi-scale and multi-directional Gabor feature vectors A(·) are used to describe each voxel to capture image texture information to 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, bandwidth 0.05, and the parameters of the low-frequency features are: wavelength 10-20, frequency 0.05-0.1, bandwidth 0.

4. For multi-direction, the direction range of the Gabor filter is set to be transformed from 0, π / 4, π / 2, 3π / 4 to π, so as to obtain image edge information from vertical to horizontal: In formula (3), η -1 (·) represents the mapping of node γ∈Γ to each common point n∈Ω t The inverse mapping function of .

7. The non-rigid registration method for intraoperative ultrasound and preoperative images based on topological high-order MRF according to claim 5, characterized in that: In step 3, in the smoothing term, 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 imposes a smooth constraint on the deformation field: P γη (L γ ,L η )=min(λ‖l γ -L η ‖2,c) (4) The minimum value of the binary cluster potential function is defined as the bi-norm of the vector difference of the binary cluster. In formula (4), l γ ,l η They represent the labels of two adjacent nodes γ and η in the target image T, λ = 2 represents the smoothness weight, and c = 0.1 is used to control the maximum cost.

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 3, in the topology term, the values of the angular Jacobian determinants corresponding to the four vertices in the four-element clique C4 in the two-dimensional image 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 four-element clique to the eight-element clique constraint, and the topological constraint condition is expanded from 4 dimensions to 64 dimensions; In step 3, in order to avoid the high cost of topological terms composed of high-order groups, a logarithmic function is proposed to penalize large deformations; for three-dimensional images, the topological term can be written as: In formula (5), w is the cost value when the topology is not preserved, satisfying log(J(l γ ,l η ,l μ ,l ψ )+1)<<w<1.

9. 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 is specifically: using the high-order variable terms and their corresponding assignments that cannot be the global minimum as the added value of the MRF potential function; the added value can reduce the order of the potential function and does not change the value of the variable corresponding to the minimum value of the potential function; the added value is expressed as a polynomial: In formula (6), C h represents a high-order group, k C is a real coefficient, and an exhaustive search method is used to find the high-order variable term ω for the MRF potential function.

10. 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 is used to achieve MRF non-rigid registration; the input image is processed at the ratio of 1 / 4, 1 / 2 and 1 / 1, and a three-level processing is used to achieve sparse sampling of coarse label distribution in a large search space to dense sampling of fine label distribution in a small search space; Low-level registration obtains the deformation field based on the energy function; mid-level and high-level registration obtain the deformation field through the energy function and the previous deformation, and upsample the multi-saliency map calculated at the previous resolution; The sampling strategy of multi-resolution processing transitions from sparse sampling in a wide search area to dense sampling in a limited area, thereby optimizing the registration accuracy and computational efficiency.

Images