A 3D CT / PET image registration method based on mutual information and L-BFGS optimization

Abstract

The present invention discloses a three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization, which belongs to the technical field of calculation, inference or counting. The registration method of the present invention is a two-stage cascade process of rigid coarse registration and B-spline elastic fine registration. It utilizes the precision advantage of mutual information measurement and noise reduction preprocessing to solve the problem of lack of precision of traditional multimodal image registration methods due to data modality and grayscale differences. It also adopts L-BFGS to improve the optimization algorithm, which greatly reduces the number of iterative calculations in traditional registration, gives full play to the speed advantage of the second-order convergence optimization algorithm, and greatly reduces the registration time. In summary, the registration method of the present invention shows high precision and speed advantages in the three-dimensional CT / PET multimodal image registration problem.

Description

Technical Field

[0001] The invention discloses a three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization, relates to computer image processing technology, and belongs to the technical field of calculation, estimation or counting. Background Art

[0002] Image registration technology aims to achieve optimal alignment of geometric structure, spatial position, and feature points for images of the same object or similar objects acquired under different conditions. It is widely used in important research areas such as computer vision and medical image processing. Modern medical imaging can be divided into two categories: anatomical imaging and functional imaging. Computed tomography (CT) is an anatomical imaging technique characterized by high resolution and the ability to provide anatomical morphological information of organs and tissues; positron emission computed tomography (PET) is a functional imaging technique characterized by the ability to provide functional metabolic information of organs and tissues, but with lower imaging resolution. In medical clinics, it is often necessary to register and integrate regional images of two modalities that present complementary information, CT imaging and PET imaging, to improve diagnostic capabilities and meet key needs such as lesion monitoring, treatment verification, intraoperative navigation, and target positioning.

[0003] Multimodal medical image registration is a very challenging task in the field of medical image processing. The large differences between CT images and PET images in imaging principles, grayscale distribution, resolution and edge features bring obvious difficulty to the registration task. In the traditional process of multimodal medical image registration algorithm, it mainly includes four modules: geometric transformation, image interpolation, similarity measurement, and search optimization. The geometric transformation parameters are iteratively optimized to maximize the similarity measurement, thereby obtaining the final registration result. Among them, the types of geometric transformation and image interpolation are usually determined in combination with the actual needs of clinical medicine. Geometric transformation mainly includes linear transformation and nonlinear transformation. Rigid transformation and B-spline elastic transformation are two types of transformation models with good registration effect and representativeness in these two types of transformations. Their principles and formulas are expressed as follows:

[0004] (1) 3D-rigid transformation

[0005] The relative position and distance between any two points in the physical space of the image remain constant before and after the 3D-rigid transformation, which can be decomposed into two transformations: translation and rotation.

[0006] Y=AX+b

[0007] Taking the 3D volume data of the present invention as an example, in the above formula, Y = (y1, y2, y3) and X = (x1, x2, x3) are the physical space coordinates of the voxel, A is a 3×3 rotation parameter matrix, b is a 3×1 translation parameter vector, and the matrix satisfies the following constraints, where I is the identity matrix and the transpose of matrix A is A T .

[0008] A T A=I,detA=1

[0009] If homogeneous coordinates are used, it can be written as follows:

[0010]

[0011] In this formula, r is a 3×3 rotation matrix, t is a translation transformation vector, p is a projection transformation vector, w produces a scaling transformation, and in a rigid body transformation, p = (0, 0, 0), w = 1, and the angle α of the image rotation around the i axis is i Use r (i) express.

[0012] r=r (1) r (2) r (3)

[0013]

[0014] In the above formula, the rotation angles of the three coordinate axes are represented by α1, α2, and α3 respectively.

[0015] (2) 3D-B spline elastic transformation

[0016] Define a set of control points φ in the original 3D image space i,j,k The size of the composition is n x ×n y ×n z For a three-dimensional grid, the voxel spacing of each control point in the three-dimensional direction is unified as δ. The 3D-B spline elastic transformation is divided into the following two steps. The first step is to calculate the floating-point grid coordinates (x, y, z), integer grid coordinates (i, j, k), and the fractional part of the floating-point grid coordinates (u, v, w) of the pixel point; the second step is to calculate the coordinate offset and interpolate.

[0017] The various parameters in the first step are mainly calculated using the row and column characteristics of the original three-dimensional image and the selected grid. The second step is to calculate the coordinate offset T shift See the following formula,

[0018]

[0019] In the above formula, B a (u), B b (v) B c (w) is the weight coefficient of the control point on the three components of the decimal part of the floating-point grid coordinate, B l is the lth basis function of B-spline, and its specific parameters are as follows,

[0020]

[0021] After obtaining the coordinate offset, the corresponding coordinates after deformation can be calculated as follows:

[0022] T(x′,y′,z′)=(x,y,z)+T shift (x,y,z)

[0023] T(x′, y′, z′) represents the coordinates of the transformed image mapped to the original image, (x, y, z) is the coordinates of the transformed image, and (x′, y′, z′) is the coordinates of the corresponding point on the original image. Because (x′, y′, z′) is a floating-point coordinate, the trilinear interpolation algorithm is also required to calculate its voxel value.

[0024] In the traditional multimodal registration algorithm process, the key to affecting the registration performance and effect lies in the latter two modules:

[0025] (3) Similarity measure calculation

[0026] The similarity measurement module primarily measures the accuracy of the current registration by calculating the similarity objective function between the reference image and the floating image during the registration process. The mutual information metric is recognized as one of the most reliable multimodal similarity measures. Its calculation principle is as follows: for two given images X and Y, H(X) and H(Y) are the entropies of images X and Y, respectively; H(X,Y) is the joint entropy of the two. The probability distribution of grayscale values ​​in images X and Y is represented by p(x) and p(y), respectively, and p(x,y) is defined as the joint probability distribution of the grayscale values ​​of images X and Y. Therefore, when H(X) = H(Y) = H(X,Y), the mutual information MI(X,Y) is maximized:

[0027] MI(X,Y)=H(X)+H(Y|X)=H(X)+H(Y)-H(X,Y)

[0028]

[0029] As one of the most widely used classic indicators, the mutual information similarity measurement index has the advantage of high precision, but its calculation time complexity is high, which is a key factor limiting the registration speed in multi-round iterative optimization calculations.

[0030] (4) Optimization and search algorithms

[0031] In the search optimization module, the input similarity measure objective function is optimized using a search strategy within the geometric transformation parameter space. After multiple rounds of iteration, the final geometric transformation parameters are obtained. The choice of optimization and search algorithm directly impacts the overall registration efficiency and accuracy of the model. Classic optimization search algorithms for medical image registration include the Powell search method, gradient descent method, and Newton method. The Powell search algorithm is widely used due to its strong local search capabilities and the lack of gradient calculation. However, it has a slow convergence rate and weak global search capabilities, making it prone to local extrema, which can affect the accuracy of the registration results. The gradient descent method, as a first-order optimization iterative algorithm, has a convergence rate between the Powell search method and the Newton method. However, near the optimal value of the objective function, the convergence rate decreases significantly and oscillations occur, affecting registration speed and accuracy. The Newton method exhibits second-order convergence in the registration optimization iterative calculation, significantly faster than the previous two algorithms and achieving higher registration accuracy. However, each iteration requires solving the inverse of the Hessian matrix of the objective function, resulting in high computational complexity and memory consumption.

[0032] Deep learning methods in the field of artificial intelligence are highly data-driven. However, for registration tasks, it is often difficult to collect large amounts of multimodal homologous data, and labeling requires a lot of manpower. In addition, the images to be registered usually contain multiple types of organs and tissues, making it difficult to determine the categories when calibrating image labels.

[0033] Given the computational complexity of similarity measures, the shortcomings of optimization and search methods in terms of convergence, registration accuracy, and complexity, and the limited availability of multimodal homologous data, the registration of 3D CT / PET images in different modalities can be improved from multiple perspectives, including data processing and iterative algorithm optimization, based on traditional registration algorithms. This application aims to propose a 3D CT / PET image registration method based on mutual information and L-BFGS optimization to overcome these shortcomings. Summary of the Invention

[0034] The purpose of the present invention is to address the shortcomings of the above-mentioned background technology and provide a three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization, which improves the accuracy and speed of registration, realizes the invention purpose of quickly and accurately registering three-dimensional CT / PET multimodal images, and solves the technical problems that traditional three-dimensional CT / PET multimodal image registration methods lack accuracy due to data modality and grayscale differences and the traditional registration optimization iterative calculation process is time-consuming.

[0035] The present invention adopts the following technical solutions to achieve the above-mentioned purpose:

[0036] A high-precision and fast three-dimensional CT / PET image registration method based on grayscale mutual information and L-BFGS optimization search is used in medical clinical practice to automatically register images of two modalities, anatomical imaging CT and functional imaging PET, which present complementary information. This three-dimensional CT / PET image registration method includes the following nine steps.

[0037] Step 1: Take one example of PET and one example of CT 3D volume data to be registered. Select one example of PET and CT 3D volume data as the fixed image, i.e., the 3D reference image f(X); and the other example of PET and CT 3D volume data as the image to be registered, i.e., the 3D floating image m(Y).

[0038] Step 2: Use a median filter-like algorithm to perform background noise reduction preprocessing on the 3D volume data of PET and CT modalities respectively.

[0039] Step 3: Using the 3D volume data of the 3D reference image f(X) as a reference, perform the following initialization preprocessing on the 3D floating image m(Y):

[0040] Step 3-1, interpolating or resampling the 3D floating image m(Y) so that the voxel spacing of the 3D floating image m(Y) is consistent with that of the 3D reference image f(X) in each dimension;

[0041] In step 3-2, the origin and direction parameters of the 3D floating image m(Y) are reset to make the origin and direction parameters of the 3D floating image m(Y) consistent with those of the 3D reference image f(X).

[0042] Step 4: The 3D reference image f(X) and the preprocessed 3D floating image m(Y) are passed into the registration model. The 3D floating image m(Y) and the 3D reference image f(X) are uniformly mapped to the same anatomical space, keeping the 3D reference image f(X) unchanged. The 3D floating image m(Y) coordinates are mapped to the 3D reference image f(X). That is, m(Y) is converted to m(X) through coordinate mapping. At this time, the three-dimensional voxel matrices of f(X) and m(X) are in the same three-dimensional RAS coordinate system. The spatial positions of the two images in the RAS coordinate system are initialized based on the three-dimensional medical image volume data parameters.

[0043] Step 5: In the geometric transformation module, the geometric transformation model T is optimized according to the L-BFGS optimizer. k The setting applies the corresponding transformation to the 3D floating image m(Y) after initialization in the RAS coordinate space m(X), and the voxel point coordinates of m(X) are converted from x to T k (x), and get the updated coordinates in the RAS coordinate system, where k is defined as the order of the transformation control parameter group, and the geometric transformation model T kIt is composed of multi-dimensional deformation parameters and adopts a staged combination of rigid transformation and B-spline elastic transformation. The geometric transformation model set by the L-BFGS optimizer includes a 3D-rigid transformation model and a B-spline elastic transformation model. Specifically, the deformation parameters of the 3D-rigid transformation model are the rotation parameter matrix A and the translation parameter vector b, and the deformation parameters of the B-spline elastic transformation model are the weight coefficients B of the control points. a (u), B b (v) B c (w) composition, when the registration is initialized, T k Set all parameters to 0.

[0044] Step 6: The spatial region where the 3D reference image f(X) is located in the RAS coordinate system is the ROI region, and the trilinear interpolation algorithm is used in the interpolation module to calculate the voxel grayscale value m(T k (X)).

[0045] Step 7: Calculate the grayscale image m(T) of the 3D reference image f(X) and the 3D floating image after geometric transformation in the ROI area in the similarity measurement module. k (X)) and take the negative mutual information index MI as the objective function S(f(X),m(T k (X))). In the registration initialization stage, the mutual information index MI between the 3D reference image f(X) and the 3D floating image after geometric transformation is calculated and the negative is taken as the objective function S(f(X),m(X)).

[0046] Step 8: In the search optimization module, L-BFGS optimization algorithm is used to measure the mutual information similarity objective function S(f(X),m(X)) or S(f(X),m(T k (X))) performs optimization calculations, and through iterative steps 5 to 7, the objective function value converges to obtain the optimal value, and the final geometric transformation parameter matrix T is output. n , wherein, according to the registration iterative process, the geometric transformation model adopts a two-stage transformation model of 3D-rigid transformation and B-spline elastic transformation respectively. In the 3D-rigid transformation registration stage, the 3D-rigid transformation deformation parameter matrix is ​​optimized using the L-BFGS optimization algorithm. When the L-BFGS optimization algorithm converges, the B-spline elastic transformation registration stage is entered, and the B-spline elastic transformation deformation parameter matrix is ​​further optimized using the L-BFGS optimization algorithm. The L-BFGS optimization algorithm used in the present invention is expressed as:

[0047] The Limited-memory BFGS (L-BFGS) algorithm is a variant of the BFGS algorithm and belongs to the class of quasi-Newton optimization algorithms. This invention aims to address the high memory cost required to store and update the approximate estimate H of the inverse Hessian matrix of the objective function during the optimization process of determining the optimal spatial transformation parameters using rigid and elastic two-stage registration. Furthermore, when using the mutual information metric MI, which has the advantage of precision in measuring the similarity of multimodal medical images, as the objective function, its computational complexity must be considered. In particular, when there are many transformation parameters to be determined in the elastic registration stage, the L-BFGS algorithm can effectively improve the convergence speed compared to other existing iterative optimization methods, leveraging the metric's precision while reducing the time and memory cost of iterative calculations. In addition, unlike other existing technologies that improve the memory cost of the traditional L-BFGS optimization algorithm, the present invention adopts different parameter boundary constraints for the rigid and elastic stages of registration implementation. The rigid registration stage mainly constrains the translation variables of the three degrees of freedom, namely the translation parameter vector b, to be controlled within the image boundary; the elastic registration stage mainly constrains the weight coefficients of the control points in the B-spline transformation to be within the specified interval (0, θ), where θ is freely controllable and is selected according to the parameter range of the coordinate offset threshold setting, so that the image will not be excessively distorted locally. In this algorithm, the approximate estimate H of the inverse of the Hessian matrix is ​​only based on the parameter space of the objective function after MI is negated and the last m steps of the gradient space, and the geometric transformation model T in the iteration k Corresponding to the parameters of the Trans function, the parameters of the Trans function are represented as (x1, x2, ..., x n ), the Trans function parameter vector composed of the Trans function parameters is recorded as x=(x1,x2,...,x n ), calculate the gradient vector of the objective function under the current deformation parameters based on the Trans function parameter vector, and the gradient vector is recorded as Define the parameter vector along the optimization path single step update amount as s k =x k+1 -x k , the gradient vector update amount is y k =g k+1 -g k , the parameter update vector s and gradient update vector y obtained in the last m iterations of the iterative process are retained in the ring buffer to reduce the space and time cost of each iterative calculation. The principle and formula are as follows:

[0048] Using the inverse formula of the BFGS algorithm, we can get the iterative formula for the approximate estimate H of the inverse of the Hessian matrix of the objective function, which is: in, I is the identity matrix. In each iteration, it mainly saves y under the conditions of k=1,2,...,k-1 k ,s k , combined with the initial H0, we can get H k ; Then perform the limited memory approximation, only take the latest m parameter update vectors s and the latest m gradient update vectors y for approximate derivation, which is the limited memory obtained by the L-BFGS algorithm for H k Approximation.

[0049] In the present invention, the process of calculating the optimal geometric transformation parameters by the L-BFGS algorithm is as follows:

[0050] ①Select an x0,m, β′<β<1 and a symmetric positive definite starting matrix, k is set to 0;

[0051] ②Calculate d k =-H k g k , x k+1 =x k +α k d k , where the step size α k Satisfies the Wolfe criterion: Initial setting step size α k =1;

[0052] ③ Order use The value pair will update H0 times, among which H k The calculation formula is:

[0053]

[0054] ④ Let k:=k+1 and return to step ② to iterate;

[0055] ⑤The number of iterations reaches the set upper limit or ||g k+1 ||≤ε, where ε is the minimum descent rate set in advance, that is, the target accuracy requirement has been achieved, the iteration ends, and the optimal solution x is returned n , which corresponds to the registration geometric transformation, that is, the final geometric transformation deformation parameter matrix T is obtained n .

[0056] Step 9: Apply the composite transformation parameter matrix composed of the 3D-rigid transformation deformation parameter matrix and the B-spline elastic transformation deformation parameter matrix to the 3D floating image m(Y) after spatial initialization under the RAS coordinates to obtain the registration result, and output the two-stage registration transformation displacement field of rigid registration and B-spline elastic registration and the three-dimensional volume data of the final registration result.

[0057] The present invention adopts the above technical solution and has the following beneficial effects:

[0058] (1) The present invention fully utilizes the precision advantage of mutual information similarity measurement in multimodal medical image registration and combines it with data denoising preprocessing to ensure that the algorithm has high registration accuracy when the data modality and grayscale differences are large.

[0059] (2) The registration process of the present invention adopts a two-stage geometric transformation model of rigid and B-spline elastic transformation to realize coarse and fine registration cascade, which not only meets the actual registration requirements but also reduces the computational time cost and parameter uncertainty and improves the registration accuracy.

[0060] (3) The improved L-BFGS optimization algorithm adopted in the present invention greatly reduces the number of iterative calculations in traditional registration, fully demonstrates the second-order convergence registration speed advantage, and further significantly reduces the registration time consumption.

[0061] In summary, the present invention adopts the above technical solutions to greatly improve the accuracy and speed of three-dimensional CT / PET multimodal image registration in terms of algorithm performance. When data is limited, it improves the accuracy and robustness advantages of traditional registration algorithms while achieving a speed second only to typical learning algorithms. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 Schematic diagram of the flow of the three-dimensional CT / PET image registration algorithm of the present invention.

[0063] Figure 2(a1) to Figure 2(e3) Figure 2 is a two-stage registration result diagram of the rigid transformation and B-spline elastic transformation of the present invention, wherein Figure 2(a1), Figure 2(a2) and Figure 2(a3) are the cross-sectional, coronal and sagittal images of the three-dimensional displacement field result of the first stage rigid registration, respectively; Figure 2(b1), Figure 2(b2) and Figure 2(b3) are the cross-sectional, coronal and sagittal images of the three-dimensional displacement field of the second stage B-spline elastic registration, respectively; Figure 2(c1), Figure 2(c2) and Figure 2(c3) Figures 2(d1), 2(d2) and 2(d3) are the two-dimensional cross-sectional images of the original three-dimensional data of the CT modality, namely the cross-sectional, coronal and sagittal planes; Figures 2(e1), 2(e2) and 2(e3) are the two-dimensional cross-sectional images of the final results of the two-stage registration of PET data after CT / PET dual-modality registration, namely the cross-sectional, coronal and sagittal planes.

[0064] Figure 3(a) to Figure 3(d)Figure 3(a) and (b) are three-dimensional visualizations of the spatial position relationship of the original CT / PET volume data; Figure 3(c) and (d) are three-dimensional visualizations of the spatial position relationship of the registered CT / PET volume data.

[0065] Figure 4 This is a time consumption comparison chart of the L-BFGS algorithm with improved search range constraints used in the iterative calculation of the two-stage registration of CT / PET dual-modality volume data in the present invention and other existing classical optimization algorithms. DETAILED DESCRIPTION

[0066] The present invention is further illustrated below with reference to specific examples. It should be understood that these examples are only used to illustrate the present invention and are not used to limit the present invention. After reading the present invention, modifications of various equivalent forms of the present invention made by those skilled in the art all fall within the scope defined by the claims attached to this application.

[0067] Example: The overall process of a high-precision and fast registration method for three-dimensional CT / PET images based on grayscale mutual information and L-BFGS optimization search is as follows: Figure 1 As shown, the specific steps include:

[0068] Step 1: Take one example of the 3D volume data to be registered, one from each of the PET and CT modalities. Select the CT example as the fixed image f(X), i.e., the 3D reference image f(X); and the PET example as the 3D floating image m(Y), i.e., the image to be registered.

[0069] Step 2: Perform background noise reduction preprocessing on the 3D volume data of PET and CT modalities respectively.

[0070] In step 3, the 3D volume data of the fixed image is used as a reference. CT images usually have a higher resolution, so the PET floating image needs to be interpolated to make the voxel spacing consistent with the fixed image. In addition, the origin and direction parameters also need to be adjusted to be consistent with the fixed image.

[0071] In step 4, the fixed image f(X) and the preprocessed 3D floating image m(Y) are passed into the registration model. The 3D floating image m(Y) and the 3D reference image f(X) are uniformly mapped to the same anatomical space X. f(X) remains unchanged, and m(Y) is converted to m(X). The 3D geometric center points of the two images are aligned based on the volume data parameters to complete the spatial position initialization in the RAS coordinate system.

[0072] Step 5: In the geometric transformation module, the geometric transformation model T is optimized according to the L-BFGS algorithm. k Setting (registration initialization T kEach parameter is set to 0), and the corresponding geometric transformation is applied to the 3D floating image m(X). The voxel coordinates of m(X) are converted from x to T k (x), and get the updated coordinates in the RAS coordinate system, where k is defined as the order of the transformation control parameter group, and the geometric transformation T k It is composed of multi-dimensional deformation parameters and adopts a staged combination of rigid transformation and B-spline elastic transformation. The geometric transformation model set by the L-BFGS optimizer includes a 3D-rigid transformation model and a B-spline elastic transformation model. Specifically, the deformation parameters of the 3D-rigid transformation model are the rotation parameter matrix A and the translation parameter vector b, and the deformation parameters of the B-spline elastic transformation model are the weight coefficients B of the control points. a (u), B b (v) B c (w) composition, when the registration is initialized, T k Set all parameters to 0.

[0073] Step 6: In the interpolation module, the trilinear interpolation algorithm is used to calculate the voxel grayscale value m(T k (X)), the spatial region where the 3D reference image f(X) is located in the RAS coordinate system is the ROI region.

[0074] Step 7: Calculate the grayscale image m(T) of the fixed image f(X) and the 3D floating image after geometric transformation in the ROI area in the similarity measurement module. k (X)) and take the negative mutual information index MI as the objective function S(f(X),m(T k (X))). In the registration initialization stage, the mutual information index MI between the 3D reference image f(X) and the 3D floating image after geometric transformation is calculated and the negative is taken as the objective function S(f(X),m(X)).

[0075] Step 8: In the search optimization module, L-BFGS optimization algorithm is used to measure the mutual information similarity objective function S(f(X),m(X)) or S(f(X),m(T k (X))) performs optimization calculations, and through iterative steps 5 to 7, the objective function value converges to obtain the optimal value, and the final geometric transformation parameter matrix T is output. n , wherein, according to the registration iterative process, the geometric transformation model adopts a two-stage transformation model of 3D-rigid transformation and B-spline elastic transformation respectively. In the 3D-rigid transformation registration stage, the 3D-rigid transformation deformation parameter matrix is ​​optimized using the L-BFGS optimization algorithm. When the L-BFGS optimization algorithm converges, the B-spline elastic transformation registration stage is entered, and the B-spline elastic transformation deformation parameter matrix is ​​further optimized using the L-BFGS optimization algorithm. The L-BFGS optimization algorithm used in the present invention is expressed as:

[0076] The memory-limited BFGS (L-BFGS) algorithm is a variant of the BFGS algorithm and belongs to the class of quasi-Newton optimization algorithms. The present invention aims to address the high memory cost required to store and update the approximate estimate H of the inverse Hessian matrix of the objective function during the optimization process of determining the optimal spatial transformation parameters using rigid and elastic two-stage registration. Furthermore, when using the mutual information metric MI, which has the advantage of precision in measuring the similarity of multimodal medical images, as the objective function, its computational complexity must be considered. In particular, when there are many transformation parameters to be determined in the elastic registration stage, the L-BFGS algorithm employed can effectively improve convergence speed compared to other existing iterative optimization methods, leveraging the metric's precision while reducing the time and memory cost of iterative calculations. In addition, unlike other existing technologies that improve the memory cost of the traditional L-BFGS optimization algorithm, the present invention adopts different parameter boundary constraints for the rigid and elastic stages of registration implementation. The rigid registration stage mainly constrains the translation variables of the three degrees of freedom to be controlled within the image boundary; the elastic registration stage mainly constrains the weight coefficients in the B-spline transformation within the specified interval (0, θ) based on the parameter range of coordinate offset and other parameters, so that the local image will not be excessively distorted. In this algorithm, the approximate estimate H of the inverse of the Hessian matrix is ​​only based on the parameter space of the objective function after MI is negated and the last m steps of the gradient space. The parameters of the geometric transformation Trans function model in the iteration are represented by (x1, x2, ..., x n ), the parameter vector is recorded as x=(x1,x2,...,x n ), the gradient vector is recorded as Define the parameter vector along the optimization path single step update amount as s k =x k+1 -x k , the gradient vector update amount is y k =g k+1 -g k , retain the last m parameter update vectors s and gradient update vectors y in the ring buffer to reduce the space and time cost of each iterative calculation. The principle and formula are as follows:

[0077] Using the inverse formula of the BFGS algorithm, we can get the iterative formula for the approximate estimate H of the inverse of the Hessian matrix of the objective function S(f(X),m(X)), which is: in I is the identity matrix. In each iteration, it mainly saves y under the conditions of k=1,2,...,k-1 k ,s k , combined with the initial H0, we can get H k; Then perform the limited memory approximation, and only take the nearly m parameter update vectors s and gradient update vectors y for approximate derivation, which is the limited memory obtained by the L-BFGS algorithm for H k Approximation.

[0078] In the present invention, the process of calculating the optimal geometric transformation parameters by the L-BFGS algorithm is as follows:

[0079] ①Select an x0,m, β′<β<1 and a symmetric positive definite starting matrix. k is set to 0.

[0080] ②Calculate d k =-H k g k , x k+1 =x k +α k d k , where the step size α k Satisfies the Wolfe criterion: Initial setting step size α k =1.

[0081] ③ Order use The value pair will update H0 times. Among them, H k The calculation formula is:

[0082]

[0083] ④ Let k:=k+1 and return to step ② to iterate.

[0084] ⑤The number of iterations reaches the set upper limit or ||g k+1 ||≤ε, where ε is the minimum descent rate set in advance, that is, the target accuracy requirement has been achieved, the iteration ends, and the optimal solution x is returned n , which corresponds to the registration geometric transformation, that is, the final transformation parameter matrix T is obtained by conversion n .

[0085] Among them, H k+1 is the k-1th iteration value of the inverse Hessian matrix of the objective function, H0 is the initial value of the inverse Hessian matrix of the objective function, V k Transpose the matrix, for The transposed matrix of for The transposed matrix of for The transposed matrix of s k 、 is the kth and the first step of the geometric transformation model deformation parameter matrix along the optimization path Step 1 Step update amount, s k =x k+1 -x k , x k 、x k+1 is the kth and k+1th iteration value of the geometric transformation model deformation parameter matrix, y k is the kth iteration value of the gradient vector update, y k =g k+1 -g k , g k+1 、g k is the k+1th and kth iteration value of the gradient vector.

[0086] Step 9: Combine the final rigid transformation and B-spline elastic transformation two-stage combined parameter matrix T n It is applied to the 3D floating image m(Y) after spatial initialization under RAS coordinates to obtain the registration result, and outputs the two-stage registration transformation displacement field of rigid registration and B-spline elastic registration and the final registration result 3D volume data.

[0087] Effect evaluation

[0088] In order to verify the effectiveness of the three-dimensional CT / PET image registration algorithm disclosed in the present invention, a dataset of image pairs of two modalities of human three-dimensional PET-CT is used to verify the significant advantages of the method disclosed in the present invention in algorithm performance such as registration accuracy and speed. Figure 2(a1) to Figure 2(e3) This is the result of the two-stage registration of the rigid transformation and B-spline elastic transformation of the present invention. Figure 2(a1) to Figure 2(a3) They are the cross-sectional, coronal and sagittal images of the three-dimensional displacement field results of the first stage of rigid registration, Figure 2(b1) to Figure 2(b3) They are the transverse, coronal and sagittal images of the second-stage B-spline elastic registration of the 3D displacement field, Figure 2(c1) to Figure 2(c3) They are the two-dimensional cross-sectional images of the original three-dimensional data of the CT modality, the coronal plane and the sagittal plane, Figure 2(d1) to Figure 2(d3) They are the two-dimensional cross-sectional images of the original three-dimensional data of the PET modality, the coronal plane and the sagittal plane, Figure 2(e1) to Figure 2(e3) The following are the two-dimensional cross-sectional images of the final results of the two-stage registration of PET data after CT / PET dual-modality registration: transverse, coronal, and sagittal planes. It can be seen that the rigid and elastic registrations are performed in stages, effectively reducing the coupling and parameter dimensionality of the two types of geometric transformation parameters, thereby effectively shortening the registration time. Furthermore, the second-order convergence advantage of the L-BFGS optimization algorithm in the present invention further reduces the iterative calculation time of the single-stage registration. Figure 3(a) to Figure 3(d)The positional relationship of the three-dimensional CT-PET volume data before and after registration in the RAS unified coordinate space clearly demonstrates the high registration accuracy of the method of the present invention. Figures 3(a) and 3(b) are the three-dimensional visualization results of the spatial positional relationship of the original CT / PET volume data; Figures 3(c) and 3(d) are the three-dimensional visualization results of the spatial positional relationship of the CT / PET volume data after registration. Figure 4 This comparison of the time required for the improved L-BFGS optimization algorithm used in the two-stage iterative calculations of the present invention's registration process with other existing classical optimization algorithms demonstrates the speed and performance advantages of the present method. In summary, the present invention utilizes the above technical solutions to significantly improve the accuracy and speed of 3D CT / PET multimodal image registration.

[0089] It should be noted that the above embodiments are merely examples of specific implementations of the present invention and are not intended to limit the scope of protection of the present invention. Equivalent replacements or substitutions made on the basis of the above technical solutions all fall within the scope of protection of the present invention.

Claims

1. A three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization, characterized in that: The steps include: Step 1: Acquire one set of PET modality 3D volume data and one set of CT modality 3D volume data, use one set of 3D volume data as a 3D reference image, and use the other set of 3D volume data as a 3D floating image; Step 2: performing background noise reduction preprocessing on the three-dimensional volume data of the 3D reference image and the three-dimensional volume data of the 3D floating image respectively; Step 3, using the three-dimensional volume data of the 3D reference image after background noise reduction preprocessing as a reference, performing resolution and direction preprocessing on the 3D floating image after background noise reduction preprocessing; Step 4: Initialize the spatial positions of the 3D reference image after background noise reduction and the 3D floating image after resolution and direction preprocessing in the RAS coordinate system; Step 5, performing geometric transformation on the 3D floating image after initialization of the spatial position under the RAS coordinates according to the geometric transformation model deformation parameter matrix obtained by the L-BFGS optimization algorithm; Step 6: Taking the spatial region where the 3D reference image initialized by the spatial position under the RAS coordinates is located as the region of interest, obtaining a grayscale image of the 3D floating image after the geometric transformation in step 5 within the region of interest; Step 7, calculating the mutual information index of the grayscale image of the 3D reference image initialized by the spatial position under the RAS coordinates and the 3D floating image after the geometric transformation in step 5 within the region of interest; Step 8: Taking the negative mutual information index of the grayscale image of the 3D reference image initialized by the spatial position under the RAS coordinates and the 3D floating image after the geometric transformation in step 5 within the region of interest as the objective function, repeat steps 5 to 7, using the L-BFGS optimization algorithm in the 3D-rigid transformation registration stage to obtain the 3D-rigid transformation deformation parameter matrix when the objective function converges, and using the L-BFGS optimization algorithm in the B-spline elastic transformation registration stage to obtain the B-spline elastic transformation deformation parameter matrix when the objective function converges; In step 9, the 3D-rigid transformation deformation parameter matrix when the objective function converges and the B-spline elastic transformation deformation parameter matrix when the objective function converges are combined into the geometric transformation model deformation parameter matrix obtained by the L-BFGS optimization algorithm, and the registration displacement fields of the 3D-rigid registration stage and the B-spline elastic registration stage and the final registration results are output.

2. The three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization according to claim 1, characterized in that: The specific method of performing resolution and direction preprocessing on the 3D floating image after background noise reduction preprocessing in step 3 is: Step 3-1, interpolating or resampling the 3D floating image after background noise reduction preprocessing until the voxel spacing of the 3D floating image is consistent with that of the 3D reference image; Step 3-2: Reset the origin parameters and direction parameters of the interpolated or resampled 3D floating image to the origin parameters and direction parameters of the 3D reference image.

3. The three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization according to claim 1, characterized in that: The specific method of step 4 of initializing the spatial positions of the 3D reference image after background noise reduction and the 3D floating image after resolution and direction preprocessing in the RAS coordinate system is as follows: keeping the 3D reference image after background noise reduction unchanged, mapping the 3D floating image after resolution and direction preprocessing to the 3D reference image after background noise reduction, and obtaining the image of the 3D floating image after spatial position initialization in the RAS coordinate system.

4. The three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization according to claim 1, characterized in that: The specific method of obtaining the grayscale image of the 3D floating image after the geometric transformation in step 5 in the region of interest in step 6 is: calculating the voxel grayscale value of the 3D floating image after the geometric transformation in step 5 in the region of interest.

5. The three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization according to claim 1, characterized in that: The 3D-rigid transformation deformation parameter matrix in step 8 includes a rotation parameter matrix and a translation parameter vector, and the B-spline elastic transformation deformation parameter matrix includes weight coefficients of each control point.

6. The three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization according to claim 5, characterized in that: The specific method of using the L-BFGS optimization algorithm in step 8 to obtain the 3D-rigid transformation deformation parameter matrix when the objective function converges is: The deformation parameter matrix of the geometric transformation model is x, x=(x1,x2,...,x n ), x1, x2, ..., x n are n deformation parameters consisting of the rotation parameter matrix of the 3D-rigid transformation, the translation parameter vector of the 3D-rigid transformation, and the weight coefficients of each control point of the B-spline elastic transformation; Let the gradient vector be g, According to the m parameters closest to the current iteration number, the vector iteration value and the gradient update vector iteration value are updated, and the expression in the form of the inverse formula of the BFGS algorithm is used. Iteratively updating the inverse of the Hessian matrix of the objective function; When the number of iterations reaches the upper limit or the gradient vector meets the accuracy requirement, the geometric transformation model deformation parameter matrix corresponding to the current iteration value of the inverse of the objective function Hessian matrix is ​​used as the 3D-rigid transformation deformation parameter matrix when the mutual information similarity measure is maximized; Among them, S(f(X),m(T k (X))) is the objective function, f(X) is the 3D reference image initialized with the spatial position under RAS coordinates, m(T k (X)) 3D floating image T after geometric transformation in the kth iteration k (X) Grayscale image within the region of interest, H k+1 is the k+1th iteration value of the inverse Hessian matrix of the objective function, H0 is the initial value of the inverse Hessian matrix of the objective function, s k is the kth iteration update of the geometric transformation model deformation parameter matrix along the optimization path, s k =z k+1 -z k , z k 、z k+1 is the kth and k+1th iteration value of the geometric transformation model deformation parameter matrix, y k is the kth iteration value of the gradient vector update, y k =g k+1 -g k , g k+1 、g k is the k+1th and kth iteration value of the gradient vector, For s k Transpose the matrix, is the deformation parameter matrix of the geometric transformation model along the optimization path sequence The iterative update amount, V k Transpose the matrix, for The transposed matrix of for The transposed matrix of for The transposed matrix of for The transposed matrix of for The transposed matrix of 7. The three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization according to claim 6, characterized in that: The specific method of using the L-BFGS optimization algorithm in step 8 to obtain the B-spline elastic transformation deformation parameter matrix when the objective function converges is the same as the specific method of using the L-BFGS optimization algorithm to obtain the 3D-rigid transformation deformation parameter matrix when the objective function converges.

8. The three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization according to claim 6, characterized in that: In the process of using the L-BFGS optimization algorithm to obtain the 3D-rigid transformation deformation parameter matrix when the objective function converges in step 8, the translation parameter vector is constrained to be within the image boundary of the 3D floating image after the spatial position is initialized under the RAS coordinate.

9. The three-dimensional CT / PET image registration method based on mutual information and L-BFGS optimization according to claim 7, characterized in that: In the process of using the L-BFGS optimization algorithm to obtain the B-spline elastic transformation deformation parameter matrix when the objective function converges in step 8, the weight coefficient of each control point is constrained to be within the specified interval (0, θ), and θ is selected according to the parameter range of the coordinate offset setting threshold during the B-spline elastic transformation process.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the method according to claim 1 is implemented.

Images