CN114757982A - Registration method and device applied to liver ablation postoperative evaluation - Google Patents

Registration method and device applied to liver ablation postoperative evaluation Download PDF

Info

Publication number
CN114757982A
CN114757982A CN202210375797.8A CN202210375797A CN114757982A CN 114757982 A CN114757982 A CN 114757982A CN 202210375797 A CN202210375797 A CN 202210375797A CN 114757982 A CN114757982 A CN 114757982A
Authority
CN
China
Prior art keywords
field
ablation
contraction
tissue
equation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210375797.8A
Other languages
Chinese (zh)
Inventor
艾丹妮
刘定坤
杨健
范敬凡
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202210375797.8A priority Critical patent/CN114757982A/en
Publication of CN114757982A publication Critical patent/CN114757982A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30056—Liver; Hepatic
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30096—Tumor; Lesion

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

The registration method and the registration device applied to liver ablation postoperative evaluation can prevent mismatching of an ablation area and a tumor, quantify external contraction of the ablation area, compensate internal contraction of the ablation area, align a microwave ablation preoperative image with a microwave ablation postoperative image, and evaluate an ablation boundary. The method comprises the following steps: (1) acquiring 3D images before and after microwave ablation as a moving image and a fixed image; performing optimization, regularization and alternating update to quantify the shrinkage field outside the ablation zone and estimate the respiratory motion field; solving the BHE to calculate a steady-state temperature field based on the ablation needle position and power, solving the TWE to obtain a steady-state contraction field based on the steady-state temperature field and the ablation zone external contraction field, merging the steady-state contraction field and the ablation zone external contraction field and obtaining a total tissue contraction field; (4) the moving image is subjected to sequential correction of respiratory motion and tissue contraction is compensated to obtain a compensated image.

Description

Registration method and device applied to liver ablation postoperative evaluation
Technical Field
The invention relates to the technical field of medical image processing, in particular to a registration method applied to post-liver ablation evaluation and a registration device applied to post-liver ablation evaluation.
Background
The conventional registration method (LogDemons) aims at finding a differential isoembryo transform,
Figure BDA0003590332050000011
it aligns the moving image (m) with the fixed image (f). The spatial transformation phi is encoded by the differential homeotic velocity field v by exponentials. Phi is exp (v). Generally, v is optimized by minimizing the energy function:
Figure BDA0003590332050000012
Figure BDA0003590332050000013
wherein similarity criteria
Figure BDA0003590332050000014
Measure f,
Figure BDA0003590332050000015
Figure BDA0003590332050000016
The similarity of (c). Rd(v) Is through a Gaussian kernel σvV is smoothed to ensure a regularization criterion for the differential homoembryology of the velocity field v. Differential homomorphism is only applicable to homologous tissue that has consistent gray scale in both images and can establish pixel level correspondence. However, in the evaluation of the ablation boundary, the tumor becomes inconsistent with the ablation zone gray due to the gray changes caused by the cytotoxic temperature, and the conventional registration method cannot correctly align the ablation zone and the tumor because they are not homologous tissues.
Disclosure of Invention
To overcome the defects of the prior art, the technical problem to be solved by the present invention is to provide a registration method applied to post-ablation-of-liver-ablation evaluation, which can prevent mismatching of an ablation region and a tumor, quantify the external contraction of the ablation region, compensate the internal contraction of the ablation region, align a pre-microwave-ablation-operation image with a post-microwave-ablation-operation image, and evaluate an ablation boundary.
The technical scheme of the invention is as follows: the registration method applied to post-liver ablation evaluation comprises the following steps:
(1) acquiring 3D images before and after microwave ablation as a moving image m and a fixed image f;
(2) performing optimization, regularization, and alternating updates to quantify shrinkage fields outside of the ablation region
Figure BDA0003590332050000021
And estimating the respiratory motion field
Figure BDA0003590332050000022
(3) Solving the Bio-heat Equation BHE (Bio-heat Equation) to calculate the steady-state temperature field T based on ablation needle position and power, and solving the thermoelastic Wave Equation TWE (thermo elastic Wave Equation) to obtain a solution based on the steady-state temperature field T and the ablation zone external contraction field
Figure BDA0003590332050000023
Steady state contraction field of
Figure BDA0003590332050000024
Merging steady state contraction fields
Figure BDA0003590332050000025
And a contraction field outside the ablation zone
Figure BDA0003590332050000026
And obtain the total tissue contraction field
Figure BDA0003590332050000027
(4) The moving image is subjected to a sequential correction of the respiratory motion and the tissue contraction is compensated to obtain a compensated image mt。
The step (2) of the invention focuses on the external ablation area, namely the area which is not reduced, the corresponding relation can be established through the gray scale of the images before and after the microwave ablation operation, the deformation is estimated and decomposed in the external ablation area, so as to quantify the external contraction and correct the respiratory motion; focusing on an internal ablation region, namely an ablation region, and finding a corresponding relation from the ablation region to an image before a microwave ablation operation through gray scale cannot be achieved, and providing a biomechanics model to compensate internal shrinkage in the internal ablation region; after continuous respiratory motion correction and tissue contraction compensation, the pre-microwave ablation and post-microwave ablation images are aligned and an ablation boundary can be evaluated; it is therefore possible to prevent mismatching of the ablation zone and the tumor, and to quantify the outer shrinkage of the ablation zone, compensate for the inner shrinkage of the ablation zone, align the pre-microwave ablation-surgery image with the post-microwave ablation-surgery image, and to evaluate the ablation boundary.
There is also provided a registration device for post-liver ablation evaluation, comprising:
an acquisition module configured to acquire pre-and post-microwave ablation 3D images as m and f;
an execution module configured to perform optimization, regularization, and alternating update to quantify an outer contraction field of an ablation zone
Figure BDA0003590332050000031
And estimating respiratory motion field
Figure BDA0003590332050000032
A solving module configured to solve the BHE to calculate a steady state temperature field T based on the ablation needle position and power, solve the TWE to obtain a solution based on T and
Figure BDA0003590332050000033
steady state contraction field of
Figure BDA0003590332050000034
Merging
Figure BDA0003590332050000035
And
Figure BDA0003590332050000036
and obtain a totalTissue contraction field
Figure BDA0003590332050000037
A correction compensation module configured to perform a sequential correction of the respiratory motion of the moving image and to compensate for tissue contraction to obtain a compensated image mt。
Drawings
Fig. 1 is a flow chart of a registration method applied to post-liver ablation evaluation according to the present invention.
Detailed Description
As shown in fig. 1, the registration method applied to post-ablation evaluation of liver comprises the following steps:
(1) acquiring 3D images before and after microwave ablation as a moving image m and a fixed image f;
(2) performing optimization, regularization, and alternating updates to quantify shrinkage fields outside of the ablation region
Figure BDA0003590332050000038
And estimating the respiratory motion field
Figure BDA0003590332050000039
(3) Solving the Bio-thermal Equation BHE (Bio-heat equalization) to calculate the steady state temperature field T based on the ablation needle position and power, solving the thermoelastic Wave Equation TWE (thermoelastic Wave equalization) ablation needle to obtain a gradient field based on the steady state temperature field T and the ablation zone outer contraction field
Figure BDA00035903320500000310
Steady state contraction field of
Figure BDA00035903320500000311
Merging steady state contraction fields
Figure BDA00035903320500000312
And a contraction field outside the ablation zone
Figure BDA00035903320500000313
And obtain the total tissue contraction field
Figure BDA00035903320500000314
(4) The moving image is subjected to a sequential correction of the respiratory motion and the tissue contraction is compensated to obtain a compensated image mt。
The step (2) of the invention focuses on the external ablation area, namely the area which is not reduced, the corresponding relation can be established through the gray scale of the images before and after the microwave ablation operation, the deformation is estimated and decomposed in the external ablation area, so as to quantify the external contraction and correct the respiratory motion; the step (3) focuses on the internal ablation region, namely the ablation region, and the corresponding relation cannot be found from the ablation region to the image before the microwave ablation operation through the gray scale, and a biomechanics model is provided for compensating the internal shrinkage in the internal ablation region; after continuous respiratory motion correction and tissue contraction compensation, the pre-microwave ablation and post-microwave ablation images are aligned and the ablation boundary can be evaluated; it is therefore possible to prevent mismatching of the ablation zone and the tumor, and to quantify the outer shrinkage of the ablation zone, compensate for the inner shrinkage of the ablation zone, align the pre-microwave ablation-surgery image with the post-microwave ablation-surgery image, and to evaluate the ablation boundary.
Preferably, in said step (2), the optimization is performed using error motion source discarding, wherein only the composite motion of respiratory motion and tissue contraction is estimated, thereby avoiding incorrect matching; defining an area omegafThis is to label the liver region ΩlAblation zone region omega of (2)cObtained given the velocity field v in the ith iterationiOnly in the region omegafIn the calculation of
Figure BDA0003590332050000041
And uses effective second order minimization to obtain the update field ui：
Figure BDA0003590332050000042
Wherein v is0Is the constant transformation, and the transformation is carried out,
Figure BDA0003590332050000043
and is the gradient at x, the maximum step size controlled by parameter/: u. ofi‖ui(x)‖≤l。
Preferably, in step (2), the regularization is a multi-source motion decomposition that quantifies tissue shrinkage by extracting shrinkage components from the composite motion; define a region omegalcWhich is the peripheral region of the ablation zone, the region omegalcIncluding respiratory motion and tissue contraction; using a Gaussian kernel σuSmooth update field uiTo ensure uiThe differential homoblast of (A) and (B); by using the HD method in the region omegalcMiddle decomposition of uiAnd split it into two sources, which are non-divergent breath-sample update fields that can estimate respiratory motion
Figure BDA0003590332050000051
And a non-rotation shrinkage sample update field that can quantify tissue shrinkage
Figure BDA0003590332050000052
ui←ui*σu,
Figure BDA0003590332050000053
Preferably, in the step (2),
Figure BDA0003590332050000054
and
Figure BDA0003590332050000055
is independently updated to obtain the respiration-like velocity field of the next iteration
Figure BDA0003590332050000056
And shrinkage velocity field of the like
Figure BDA0003590332050000057
Figure BDA0003590332050000058
Figure BDA0003590332050000059
Wherein
Figure BDA00035903320500000510
And
Figure BDA00035903320500000511
is identity transformation, BCH (·,. cndot.) is a Beck-Campbell-Hostaff formula, and can effectively calculate the updated velocity field in the logarithmic domain;
Figure BDA00035903320500000512
and
Figure BDA00035903320500000513
respectively representing respiratory motion
Figure BDA00035903320500000514
And tissue contraction
Figure BDA00035903320500000515
Figure BDA00035903320500000516
Combined with the effective BCH formula to obtain the next optimized composite velocity field vi+1：
Figure BDA00035903320500000517
Preferably, in the step (2), in the registration, the following steps are alternately performed: optimization procedureTo obtain an update field uiThe regularization step is to normalize uiDecomposition into breath-like update fields
Figure BDA00035903320500000518
And a shrinkage sample update field
Figure BDA00035903320500000519
The updating step obtains a new breath sample velocity field
Figure BDA00035903320500000520
And shrinkage sample velocity field
Figure BDA00035903320500000521
And calculating the composite velocity field v of the next iterationi+1；
When the registration converges, we can get two combinable displacements
Figure BDA00035903320500000522
And
Figure BDA00035903320500000523
representing respiratory motion and tissue contraction, respectively:
Figure BDA00035903320500000524
wherein
Figure BDA00035903320500000525
And
Figure BDA00035903320500000526
respectively estimating the displacement, the respiratory sample velocity field and the contraction sample velocity field;
multi-source deformation is broken down into separate parts, first considering the total displacement as respiratory motion, then tissue contraction:
Figure BDA0003590332050000061
mdis a moving image corrected for respiratory motion.
Preferably, in step (3), the tissue is considered homogeneous, and the steady-state temperature field T is obtained by solving BHE:
Figure BDA0003590332050000062
Figure BDA0003590332050000063
where c and ρ are the specific heat capacity and mass density of the tissue, the left side of the equation is equal to 0, considering that the temperature field reaches steady state, with its bias derivative with respect to time being 0; the right side of the equation represents different mechanisms of heat sinking and dissipation: heat transfer by internal conduction is the thermal conductivity of tissue, and the heat exchange mechanism by metabolic heat generation, electromagnetic power sink, and capillary blood perfusion is blood density, c bIs the specific heat capacity of blood, omegabIs the blood perfusion rate, TbIs the blood vessel temperature, metabolic heat production is neglected;
the distribution in equation (9) is determined by the structural characteristics of the microwave ablation surgical ablation needle and is obtained by fitting the experimentally measured data SAR with an exponential fit r in the radial direction and a cubic polynomial fit z in the axial direction:
SAR(r,z；P)＝αeβr(c0z3+c1z2+c2z1+c3) (10)
p where is a given microwave power, alpha, beta, c0、c1、c2And c3Is the parameter to be fitted;
equation (9) is a partial differential equation relating the temperature T at the edge of the ablation needleaAnd condensation edge temperature TcAs a dirichlet boundary condition, the steady-state temperature field T is obtained by solving the equation using the finite difference FD.
Preferably, in the step (3), the steady-state contraction field is calculated by solving two:
Figure BDA0003590332050000064
Figure BDA0003590332050000065
Figure BDA0003590332050000066
wherein V Poisson's ratio, V is the evaporation amount of water, w is the relative water content, λ and μ are the volumetric thermal expansion coefficient and volumetric water vaporization contraction coefficient, respectively, and G and K are the shear modulus and the bulk modulus, respectively, in consideration of the contraction field
Figure BDA0003590332050000071
A steady state is reached, the derivative of its bias with respect to time is 0, and the left side of the equation equals 0.
Preferably, in the step (3), the shear modulus and the bulk modulus are obtained as follows:
Figure BDA0003590332050000072
wherein E is Young's modulus;
the relative water content w is related to the relaxation time T2, and is calculated by the spin-spin relaxation signal decay equation (13):
Figure BDA0003590332050000073
Where TE is the echo time, M is the spin density, S is the signal intensity value in the T2 weighted image, W is used to measure the proportion of the moisture content in the different tissues as an indication of the contraction potential of the different tissues W, W is normalized;
the water evaporation amount V is calculated by the formula (14):
Figure BDA0003590332050000074
preferably, in the step (4), merging
Figure BDA0003590332050000075
And obtain the total tissue contraction field
Figure BDA0003590332050000076
Quantification and compensation of tissue contraction:
Figure BDA0003590332050000077
using the total tissue shrinkage
Figure BDA0003590332050000078
Updating equation (8) results in the final warped image:
Figure BDA0003590332050000079
Figure BDA00035903320500000710
it will be understood by those skilled in the art that all or part of the steps in the method of the above embodiments may be implemented by hardware instructions related to a program, the program may be stored in a computer-readable storage medium, and when executed, the program includes the steps of the method of the above embodiments, and the storage medium may be: ROM/RAM, magnetic disks, optical disks, memory cards, and the like. Therefore, in correspondence with the method of the invention, the invention also comprises a registration device for post-ablation evaluation of the liver, which is generally represented in the form of functional modules corresponding to the steps of the method. The device comprises an acquisition module configured to acquire an image to be segmented;
An acquisition module configured to acquire pre-and post-microwave ablation 3D images as m and f;
an execution module configured to perform optimization, regularization, and alternating update to quantify an ablation region external contraction field
Figure BDA0003590332050000081
And estimating the respiratory motion field
Figure BDA0003590332050000082
A solving module configured to solve the BHE to calculate a steady state temperature field T based on the ablation needle position and power, solve the TWE to obtain a solution based on T and
Figure BDA0003590332050000083
steady state contraction field of
Figure BDA0003590332050000084
Merging
Figure BDA0003590332050000085
And
Figure BDA0003590332050000086
and obtain the total tissue contraction field
Figure BDA0003590332050000087
A correction compensation module configured to perform a sequential correction of the respiratory motion of the moving image and to compensate for tissue contraction to obtain a compensated image mt。
The present invention will be described in more detail below.
The initially estimated velocity field consists of a combination of multiple motion sources, with respiratory motion and tissue contraction, and a mismatch between the ablation zone and the tumor, which can lead to severe distortion. The present invention proposes v-local shrinkage decomposition (LC module) to eliminate incorrect matching and decompose complex deformations to quantify the shrinkage outside the ablation zone.
The LC module consists of three steps. The first step is optimization using error motion source dropping, where onlyThe composite motion of the respiratory motion and tissue contraction is estimated, thereby avoiding incorrect matching. In particular, a region Ω is defined fThis is to label the liver region ΩlAblation zone region omega of (2)cAnd obtaining the product. Given the velocity field v in the ith iterationiOnly in the region omegafIn the calculation of
Figure BDA0003590332050000088
And using effective second order minimization to obtain the update field ui：
Figure BDA0003590332050000091
Wherein v is0Is the constant transformation, and the transformation is carried out,
Figure BDA0003590332050000092
and is the gradient at x. Maximum step size controlled by parameter/: u. ofi‖ui(x)‖≤l
Considering that the intensity of the ablation region approximates the intensity of the lesion, an incorrect match may only occur at Ωc. Optimization procedure in omegafAnd is not at omegacIn
Figure BDA0003590332050000093
The resulting update field uiSources of error motion can be excluded, thereby preventing incorrect matching of ablation regions to the tumor.
The second step of the LC module, the regularization step, is a multi-source motion decomposition that quantifies tissue shrinkage by extracting shrinkage components from the composite motion. In particular, a region Ω is definedlcWhich is the peripheral region of the ablation zone. The region omegalcMainly involving respiratory motion and tissue contraction. Helmholtz Decomposition (HD) demonstrated that any differential homeomorphic displacement field u can be uniquely decomposed into a divergence-free field and a rotation-free field: u-ud+ucIn which there is no field of divergence
Figure BDA0003590332050000094
And no rotation field
Figure BDA0003590332050000095
Figure BDA0003590332050000096
In the regularization step, a Gaussian kernel σ is useduSmooth update field uiTo ensure uiThe differential of (A) is isoembryonal. Then, in the region Ω by using the HD method lcMiddle decomposition of uiAnd divide it into two sources, i.e. non-divergence breath sample update fields, where the respiratory motion can be estimated
Figure BDA0003590332050000097
And a rotation-free shrinkage sample update field that can quantify tissue shrinkage
Figure BDA0003590332050000098
ui←ui*σu
Figure BDA0003590332050000099
The third step of the LC module is the update.
Figure BDA00035903320500000910
And
Figure BDA00035903320500000911
is independently updated to obtain the respiration-like velocity field of the next iteration
Figure BDA00035903320500000912
And shrinkage velocity field of the like
Figure BDA00035903320500000913
Figure BDA00035903320500000914
Figure BDA00035903320500000915
Wherein
Figure BDA00035903320500000916
And
Figure BDA00035903320500000917
is an identity transformation. BCH (·,. cndot.) is a Becky-Campbell-Housdov formula that efficiently computes the updated velocity field in the logarithmic domain.
Figure BDA0003590332050000101
And
Figure BDA0003590332050000102
can respectively represent breathing movement
Figure BDA0003590332050000103
And tissue contraction
Figure BDA0003590332050000104
Figure BDA0003590332050000105
They can be combined with the effective BCH formula to obtain the next optimized composite velocity field vi+1：
Figure BDA0003590332050000106
In registration, the following steps are performed alternately: optimization procedure to obtain an update field uiThe regularization step is to normalize uiDecomposition into breath-like update fields
Figure BDA0003590332050000107
And a shrinkage sample update field
Figure BDA0003590332050000108
The updating step obtains a new breath sample velocity field
Figure BDA0003590332050000109
And shrinkage sample velocity field
Figure BDA00035903320500001010
And calculating the composite velocity field v of the next iterationi+1。
When the registration converges, two combinable displacements can be obtained
Figure BDA00035903320500001011
And
Figure BDA00035903320500001012
representing respiratory motion and tissue contraction, respectively:
Figure BDA00035903320500001013
wherein
Figure BDA00035903320500001014
And
Figure BDA00035903320500001015
respectively, estimates of displacement, respiration sample velocity field, and contraction sample velocity field.
Based on the LC module, the multi-source deformation is decomposed into independent parts. First consider the total displacement as respiratory motion, followed by tissue contraction:
Figure BDA00035903320500001016
mdIs a moving image corrected for respiratory motion. Using an LC module, respiratory motion correction may be quantified
Figure BDA00035903320500001017
External contraction of posterior ablation zone
Figure BDA00035903320500001018
The LC module estimates and resolves the deformation of the outer ablation zone region
Figure BDA00035903320500001019
To quantify the shrinkage outside the ablation zone. The invention proposes biomechanical model constraints (BM modules) to compensate for the internal shrinkage of ablation zones based on external shrinkage
Figure BDA00035903320500001020
In clinical practice, the physician typically heats the tumor at a cytotoxic temperature for a sufficient time to reach a sufficient ablation boundary. Recent studies have shown that at constant power and sufficient heating time, a steady-state thermal field and a steady-state contraction field can be achieved. Based on these studies, the BHE was solved to calculate the steady state temperature field T based on the ablation needle position and power. Furthermore, the TWE is solved to obtain a Tsum-based
Figure BDA00035903320500001021
Steady state contraction field of
Figure BDA00035903320500001022
The tissue is considered homogeneous and the steady state temperature field T can be obtained by solving for BHE:
Figure BDA0003590332050000111
Figure BDA0003590332050000112
wherein c (J.kg)-1·K-1) And rho (kg. m)-3) Is the specific heat capacity and mass density of the tissue. Taking into account that the temperature field reaches a steady state, its offset with respect to timeThe derivative is 0 and the left side of the equation equals 0. The right side of the equation represents different mechanisms of heat sinking and dissipation: heat transfer by internal conduction ((k (W.m)) -1·K-1) Is tissue heat conductivity), and metabolizes to generate heat (W.m)-3) Electromagnetic power deposition SAR (W.m)-3) And capillary blood perfusion induced heat exchange mechanism (ρ)b(kg·m-3) Is the blood density, cbIs specific heat capacity (J.kg) of blood-1·K-1)，ωb(kg·m-3·s-1) Is the blood perfusion rate). T is a unit ofb(K) Is the blood vessel temperature. Metabolic thermogenesis a is ignored because its contribution is very small relative to the other terms.
(9) Is determined by structural characteristics of the microwave ablation surgical ablation needle and is obtained by fitting experimentally measured data SAR with an exponential fit (r) in the radial direction and a cubic polynomial fit (z) in the axial direction:
SAR(r,z；P)＝αeβr(c0z3+c1z2+c2z1+c3) (10)
p where is a given microwave power. Alpha, beta, c0、c1、c2And c3Are the parameters to be fitted (see table I).
Equation (9) is a partial differential equation that requires a boundary condition to obtain a unique solution. Temperature T of the edge of the ablation needleaAnd condensation edge temperature TcAs a dirichlet boundary condition. The steady state temperature field T is obtained by solving (9) the equation using Finite Differences (FD), where the parameters are shown in table I.
Then, the steady state shrinkage field is calculated by solving for TWE:
Figure BDA0003590332050000113
Figure BDA0003590332050000114
Figure BDA0003590332050000115
wherein v is the poisson's ratio; v is the evaporation amount of water, and w is the relative water content. And λ and μ are the volumetric thermal expansion coefficient and the volumetric water vaporization contraction coefficient, respectively. G and K are shear and bulk modulus, respectively. Taking into account the shrinkage field
Figure BDA0003590332050000116
A steady state is reached, the derivative of its bias with respect to time is 0, and the left side of the equation is equal to 0.
The shear and bulk moduli can be obtained as follows:
Figure BDA0003590332050000121
wherein E is the Young's modulus.
In a recent study, the relative water content w is related to the relaxation time T2, which can be approximated by the spin-spin relaxation signal decay equation:
Figure BDA0003590332050000122
where TE is the echo time; m is the spin density, S and is the signal intensity value in the T2 weighted image. W is used to measure the proportion of the moisture content in different tissues as an indication W of the contraction potential of the different tissues. Thus, w is normalized. Different tissues
Figure BDA0003590332050000123
See table I for details.
The water evaporation V is calculated by the piecewise function proposed by Yang et al.
Figure BDA0003590332050000124
Equation (11) is also a partial differential equation, which requires oneBoundary conditions to obtain a unique solution. Using ablation zone external contraction
Figure BDA0003590332050000125
As a dirichlet boundary condition. Obtaining a steady-state shrinkage field inside an ablation zone by solving (11) using the FD method
Figure BDA0003590332050000126
The detailed parameters are shown in table I.
Merging
Figure BDA0003590332050000127
And obtain the total tissue contraction field
Figure BDA0003590332050000128
And quantify and compensate for tissue shrinkage:
Figure BDA0003590332050000129
using total tissue shrinkage
Figure BDA00035903320500001210
Updating (8) the final warped image:
Figure BDA00035903320500001211
Figure BDA00035903320500001212
given an LC module and a BM module, the proposed LC-BM elastic registration method can estimate global differential homeotropic displacement fields with local contraction, align pre-and post-microwave ablation images, and compensate tissue contraction to assess ablation boundaries.
Based on open source LogDemons realization, an LC-BM elastic registration method is realized by using an Insight Toolkit, SimpleITK, Visualization Toolkit and an ArrayFire GPU matrix library. Table i detailed parameters of the proposed LC-BM method are given.
Table I.
Parameter list for LC-BM method
Figure BDA0003590332050000131
The above description is only a preferred embodiment of the present invention, and is not intended to limit the present invention in any way, and any simple modifications, equivalent variations and modifications made on the above embodiment according to the technical spirit of the present invention are within the scope of the technical solution of the present invention.

Claims (10)

1. The registration method applied to liver ablation postoperative evaluation is characterized by comprising the following steps: which comprises the following steps:
(1) acquiring 3D images before and after microwave ablation operation as a moving image m and a fixed image f;
(2) performing optimization, regularization, and alternating updates to quantify the shrinkage field outside the ablation zone
Figure FDA0003590332040000011
And estimating respiratory motion field
Figure FDA0003590332040000012
(3) Solving a biological thermal equation BHE to calculate a steady-state temperature field T based on the ablation needle position and power, and solving a thermoelastic wave equation TWE to obtain a model based on the steady-state temperature field T and an ablation zone external contraction field
Figure FDA0003590332040000013
Steady state contraction field of
Figure FDA0003590332040000014
Merging steady state contraction fields
Figure FDA0003590332040000015
And a contraction field outside the ablation zone
Figure FDA0003590332040000016
And obtain the total tissue contraction field
Figure FDA0003590332040000017
(4) The moving image is subjected to a sequential correction of the respiratory motion and the tissue contraction is compensated to obtain a compensated image mt。
2. The registration method applied to post-liver ablation assessment according to claim 1, wherein: in the step (2), an error motion source is discarded for optimization, wherein only the composite motion of respiratory motion and tissue contraction is estimated, so as to avoid incorrect matching; defining an area omegafThis is to label the liver region ΩlAblation zone region omega of (2)cObtained given the velocity field v in the ith iterationiOnly in the region omegafMiddle calculation
Figure FDA00035903320400000110
Figure FDA00035903320400000111
And uses effective second order minimization to obtain the update field ui：
Figure FDA0003590332040000018
Wherein v is0Is the constant transformation, and the transformation is carried out,
Figure FDA0003590332040000019
and is the gradient at x, the maximum step size controlled by parameter/: u. ofi‖ui(x)‖≤l。
3. The registration method for post-liver ablation assessment according to claim 2, wherein: in the step (2)Regularization is a multi-source motion decomposition that quantifies tissue shrinkage by extracting shrinkage components from the composite motion; define a region omegalcWhich is the peripheral region of the ablation zone, the region omegalcIncluding respiratory motion and tissue contraction; using a Gaussian kernel σuSmooth update field uiTo ensure u iDifferential isoembryogenesis of (3); decomposing HD method in region omega by using HelmholtzlcMiddle decomposition of uiAnd split it into two sources, which are non-divergent breath sample update fields that can estimate respiratory motion
Figure FDA0003590332040000021
And a rotation-free shrinkage sample update field that can quantify tissue shrinkage
Figure FDA0003590332040000022
ui←ui*σu,
Figure FDA0003590332040000023
4. The registration method applied to post-liver ablation assessment according to claim 3, wherein: in the step (2), the step (c),
Figure FDA0003590332040000024
and
Figure FDA0003590332040000025
is independently updated to obtain the respiration-like velocity field of the next iteration
Figure FDA0003590332040000026
And shrinkage velocity field of the like
Figure FDA0003590332040000027
Figure FDA0003590332040000028
Figure FDA0003590332040000029
Wherein
Figure FDA00035903320400000210
And
Figure FDA00035903320400000211
is identity transformation, BCH (·,. cndot.) is a Beck-Campbell-Hostaff formula, and can effectively calculate the updated velocity field in the logarithmic domain;
Figure FDA00035903320400000212
and
Figure FDA00035903320400000213
respectively representing respiratory motion
Figure FDA00035903320400000214
And tissue contraction
Figure FDA00035903320400000215
Figure FDA00035903320400000216
Combined with the effective BCH formula to obtain the next optimized composite velocity field vi+1：
Figure FDA00035903320400000217
5. According to the claimsSolving 4 the registration method applied to evaluation after liver ablation, which is characterized in that: in the step (2), in the registration, the following steps are alternately performed: optimization procedure to obtain an update field uiThe regularization step is to normalize uiDecomposition into breath-like update fields
Figure FDA00035903320400000218
And a shrinkage sample update field
Figure FDA0003590332040000031
The updating step obtains a new breath sample velocity field
Figure FDA0003590332040000032
And shrinkage sample velocity field
Figure FDA0003590332040000033
And calculating the composite velocity field v of the next iterationi+1；
When the registration converges, we can get two combinable displacements
Figure FDA0003590332040000034
And
Figure FDA0003590332040000035
respectively representing respiratory motion and tissue contraction:
Figure FDA0003590332040000036
wherein
Figure FDA0003590332040000037
And
Figure FDA0003590332040000038
respectively estimating the displacement, the respiratory sample velocity field and the contraction sample velocity field;
multi-source deformation is broken down into separate parts, first considering the total displacement as respiratory motion, then tissue contraction:
Figure FDA0003590332040000039
mdis a moving image corrected for respiratory motion.
6. The registration method for post-liver ablation assessment according to claim 5, wherein: in said step (3), the tissue is considered homogeneous and the steady-state temperature field T is obtained by solving BHE:
Figure FDA00035903320400000310
Figure FDA00035903320400000311
where c and ρ are the specific heat capacity and mass density of the tissue, the left side of the equation is equal to 0, considering that the temperature field reaches steady state, with its bias derivative with respect to time being 0; the right side of the equation represents different mechanisms of heat sinking and dissipation: heat transfer by internal conduction is the thermal conductivity of tissue, and the heat exchange mechanism by metabolic heat generation, electromagnetic power sink, and capillary blood perfusion is blood density, cbIs the specific heat capacity of blood, omegabIs the blood perfusion rate, TbIs the blood vessel temperature, metabolic heat production is neglected;
The distribution in equation (9) is determined by the structural characteristics of the microwave ablation surgical ablation needle and is obtained by fitting the experimentally measured data SAR with an exponential fit r in the radial direction and a cubic polynomial fit z in the axial direction:
SAR(r,z；P)＝αeβr(C0z3+c1z2+c2z1+c3) (10)
p where is a given microwave power, alpha, beta, c0、c1、c2And c3Is the parameter to be fitted;
equation (9) is a partial differential equation relating the temperature T at the edge of the ablation needleaAnd condensation edge temperature TcAs a dirichlet boundary condition, the steady-state temperature field T is obtained by solving the equation using the finite difference FD.
7. The registration method for post-liver ablation assessment according to claim 6, wherein: in the step (3), the steady-state contraction field is calculated by solving TWE:
Figure FDA0003590332040000041
Figure FDA0003590332040000042
Figure FDA0003590332040000043
wherein V Poisson's ratio, V is the evaporation amount of water, w is the relative water content, λ and μ are the volumetric thermal expansion coefficient and volumetric water vaporization contraction coefficient, respectively, and G and K are the shear modulus and the bulk modulus, respectively, in consideration of the contraction field
Figure FDA0003590332040000044
A steady state is reached, the derivative of its bias with respect to time is 0, and the left side of the equation equals 0.
8. The registration method for post-liver ablation assessment according to claim 7, wherein: in the step (3), the shear modulus and the bulk modulus are obtained by:
Figure FDA0003590332040000045
Wherein E is Young's modulus;
the relative water content w is related to the relaxation time T2, and is calculated by the spin-spin relaxation signal decay equation (13):
Figure FDA0003590332040000046
where TE is the echo time, M is the spin density, S is the signal intensity value in the T2 weighted image, W is used to measure the proportion of moisture content in different tissues as an indication W of the contractile potential of the different tissues, W is normalized;
the water evaporation amount V is calculated by the formula (14):
Figure FDA0003590332040000051
9. the registration method applied to post-liver ablation assessment according to claim 8, wherein: in the step (4), merging
Figure FDA0003590332040000052
And obtain the total tissue contraction field
Figure FDA0003590332040000053
Quantification and compensation of tissue contraction:
Figure FDA0003590332040000054
using the total tissue shrinkage
Figure FDA0003590332040000055
Updating equation (8) to obtain the finalThe distorted image of (2):
Figure FDA0003590332040000056
10. be applied to registration device of liver ablation postoperative aassessment which characterized in that: it includes:
an acquisition module configured to acquire pre-and post-microwave ablation 3D images as m and f;
an execution module configured to perform optimization, regularization, and alternating update to quantify an ablation region external contraction field
Figure FDA0003590332040000057
And estimating the respiratory motion field
Figure FDA0003590332040000058
A solving module configured to solve the BHE to calculate a steady state temperature field T based on the ablation needle position and power, solve the TWE to obtain a solution based on T and
Figure FDA0003590332040000059
Steady state contraction field of
Figure FDA00035903320400000510
Merging
Figure FDA00035903320400000511
And
Figure FDA00035903320400000512
and obtain the total tissue contraction field
Figure FDA00035903320400000513
A correction compensation module configured to perform a sequential correction of respiratory motion on the moving image and to compensate for tissue contraction toObtaining a compensated image mt。
CN202210375797.8A 2022-04-11 2022-04-11 Registration method and device applied to liver ablation postoperative evaluation Pending CN114757982A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210375797.8A CN114757982A (en) 2022-04-11 2022-04-11 Registration method and device applied to liver ablation postoperative evaluation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210375797.8A CN114757982A (en) 2022-04-11 2022-04-11 Registration method and device applied to liver ablation postoperative evaluation

Publications (1)

Publication Number Publication Date
CN114757982A true CN114757982A (en) 2022-07-15

Family

ID=82330012

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210375797.8A Pending CN114757982A (en) 2022-04-11 2022-04-11 Registration method and device applied to liver ablation postoperative evaluation

Country Status (1)

Country Link
CN (1) CN114757982A (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109859833A (en) * 2018-12-28 2019-06-07 北京理工大学 The appraisal procedure and device of ablative surgery therapeutic effect
CN110910406A (en) * 2019-11-20 2020-03-24 中国人民解放军总医院 Method and system for evaluating three-dimensional space curative effect after liver tumor ablation
WO2021088747A1 (en) * 2019-11-04 2021-05-14 中国人民解放军总医院 Deep-learning-based method for predicting morphological change of liver tumor after ablation
CN113378879A (en) * 2021-05-06 2021-09-10 上海美杰医疗科技有限公司 Postoperative tumor assessment method and device and computer storage medium

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109859833A (en) * 2018-12-28 2019-06-07 北京理工大学 The appraisal procedure and device of ablative surgery therapeutic effect
WO2021088747A1 (en) * 2019-11-04 2021-05-14 中国人民解放军总医院 Deep-learning-based method for predicting morphological change of liver tumor after ablation
CN110910406A (en) * 2019-11-20 2020-03-24 中国人民解放军总医院 Method and system for evaluating three-dimensional space curative effect after liver tumor ablation
CN113378879A (en) * 2021-05-06 2021-09-10 上海美杰医疗科技有限公司 Postoperative tumor assessment method and device and computer storage medium

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
金华等: "实时三维超声造影对肝癌射频消融安全边界的评估价值", 影像研究与医学应用, vol. 4, no. 20, 15 October 2020 (2020-10-15) *

Similar Documents

Publication Publication Date Title
Klages et al. Patch‐based generative adversarial neural network models for head and neck MR‐only planning
Berger et al. Marker‐free motion correction in weight‐bearing cone‐beam CT of the knee joint
EP1405093A1 (en) Mri-based temperature mapping with error compensation
Riblett et al. Data‐driven respiratory motion compensation for four‐dimensional cone‐beam computed tomography (4D‐CBCT) using groupwise deformable registration
Ren et al. Development and clinical evaluation of a three-dimensional cone-beam computed tomography estimation method using a deformation field map
JP5016518B2 (en) Alignment apparatus and program thereof
Peressutti et al. A novel Bayesian respiratory motion model to estimate and resolve uncertainty in image-guided cardiac interventions
Udupa et al. Understanding respiratory restrictions as a function of the scoliotic spinal curve in thoracic insufficiency syndrome: a 4D dynamic MR imaging study
Yang et al. Nonrigid registration of medical image based on adaptive local structure tensor and normalized mutual information
Dai et al. Deep learning‐based motion tracking using ultrasound images
CN108171737B (en) Medical image registration method and system with incompressible organ
CN104933672B (en) Method based on quick convex optimized algorithm registration three dimensional CT with ultrasonic liver image
Espinel et al. Using multiple images and contours for deformable 3d-2d registration of a preoperative ct in laparoscopic liver surgery
Gong et al. Multiple-object 2-D–3-D registration for noninvasive pose identification of fracture fragments
Mohammadi et al. A combined registration and finite element analysis method for fast estimation of intraoperative brain shift; phantom and animal model study
Li et al. Pulmonary CT image registration and warping for tracking tissue deformation during the respiratory cycle through 3D consistent image registration
CN114757982A (en) Registration method and device applied to liver ablation postoperative evaluation
Noorda et al. Subject-specific four-dimensional liver motion modeling based on registration of dynamic MRI
Morris et al. High-spatial-resolution bone densitometry with dual-energy X-ray absorptiometric region-free analysis
Chen et al. Accurate total‐body Ki parametric imaging with shortened dynamic 18F‐FDG PET scan durations via effective data processing
CN109767410B (en) Lung CT and MRI image fusion algorithm
Thomas et al. Computer-assisted contralateral side comparison of the ankle joint using flat panel technology
Zhao et al. Chest modeling and personalized surgical planning for pectus excavatum
Koolwal et al. A fast SLAM approach to freehand 3-D ultrasound reconstruction for catheter ablation guidance in the left atrium
White et al. Modeling and incorporating cardiac‐induced lung tissue motion in a breathing motion model

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination