An ablation computing method and an ablation computing system

Abstract

An ablation calculation method and an ablation calculation system, the ablation calculation method comprising: scanning a target to be measured using a gradient echo sequence containing i different echo times, obtaining a phase map corresponding to the echo time, i being a positive integer greater than or equal to 2 (S101); selecting at least two groups of phase maps corresponding to different echo times to obtain corresponding temperature difference maps (S102); obtaining a temperature map according to the temperature difference maps (S103); calculating an ablation condition according to the temperature map (S104); and the ablation calculation system comprising an ablation calculation module capable of executing the ablation calculation method.

Description

[0001] This application claims priority to Chinese Patent Application No. 202110184184.1, filed on February 8, 2021, entitled “An Ablation Calculation Method and Ablation Calculation System,” the contents of which are incorporated herein by reference in their entirety for all purposes. Technical Field

[0002] This application relates to the field of medical devices, and more specifically, to an ablation calculation method and ablation calculation system. Background Technology

[0003] Magnetic Resonance Temperature Imaging (MRTI) enables non-invasive, real-time, in vivo monitoring of the internal temperature distribution and changes of a subject. It has important applications in the ablation monitoring process of minimally invasive and non-invasive hyperthermia treatments, such as magnetic resonance interstitial hyperthermia and focused ultrasound therapy.

[0004] One current magnetic resonance imaging method is the proton resonance frequency (PRF) displacement-based thermometry. However, in practice, it has been found that this method is significantly affected by objective environmental factors such as the uniformity of the magnetic field of the magnetic resonance coil, the non-uniform distribution of tissue magnetic susceptibility, and tissue movement caused by respiration / blood flow pulsation. Common errors caused by these factors include phase unwrapping errors, errors caused by rapid changes in magnetic susceptibility, and errors caused by motion. This can easily lead to a large discrepancy between the final temperature map and the actual temperature, rendering the temperature map meaningless.

[0005] How to perform more accurate ablation calculations based on magnetic resonance data is a technical problem that still needs to be solved in this field. Summary of the Invention

[0006] To address the aforementioned technical problems, this application provides an ablation calculation method and a related ablation calculation system.

[0007] In a first aspect of the invention, an ablation calculation method is provided, the method comprising:

[0008] The target under test is scanned using a gradient echo sequence containing i different echo times to obtain a phase map corresponding to the echo times, where i is a positive integer greater than or equal to 2;

[0009] Select at least two sets of phase diagrams corresponding to different echo times to obtain the corresponding temperature difference diagrams;

[0010] A temperature map is obtained based on the temperature difference map;

[0011] The ablation status (for each pixel) is calculated based on the temperature map.

[0012] Furthermore, in this method, the temperature difference map is obtained as follows: the phase difference map at any given time is obtained by subtracting the phase map at a reference time from the phase map at any given time; at least one phase difference map corresponding to an echo time is selected as a reference phase difference map; the phase difference maps to be calibrated corresponding to other echo times are calibrated based on the reference phase difference map to obtain a calibrated phase difference map; the echo time corresponding to the reference phase difference map is less than the echo time corresponding to the calibrated phase difference map; the temperature difference map at that time is calculated using the reference phase difference map and the calibrated phase difference map.

[0013] The echo time corresponding to at least one reference phase difference map does not exceed one of the following: 18ms, 17ms, 16ms, 15ms, 14ms, 13ms, 12ms, 11ms, 10ms, 9ms, 8ms, 7ms, 6ms, 5ms or 4ms.

[0014] Furthermore, in this method, the above calibration is performed as follows:

[0015] Based on the proportional relationship between the phase difference diagram and the echo time, the estimated value of the phase difference diagram to be calibrated is calculated based on the echo time and the reference phase difference diagram.

[0016] Using the estimated value of the phase difference map to be calibrated, the phase difference map to be calibrated is unwrapped according to the phase periodicity to obtain the calibrated phase difference map.

[0017] Optionally, the ablation calculation method of the present invention further includes a step of eliminating phase drift caused by the magnetic resonance system (e.g., B0 drift error); further, in the step of eliminating phase drift caused by the magnetic resonance system, several regions with stable physical temperature and uniform tissue are selected as thermal reference points, and phase drift correction is performed by subtracting the average phase difference of the thermal reference points from each phase difference image or subtracting the average temperature change of the thermal reference points from the temperature change image.

[0018] Optionally, the ablation calculation method of the present invention further includes a step of correcting errors caused by magnetic susceptibility. The magnetic susceptibility correction step is performed on a phase difference diagram or a temperature difference diagram, and this step includes:

[0019] The steps for performing magnetic susceptibility correction on a temperature difference graph include:

[0020] A first temperature map is obtained based on the reference phase difference map, and a corresponding second temperature map is obtained based on the calibrated phase difference map.

[0021] Determine whether the absolute value of the difference between the temperature value corresponding to each pixel in the second temperature map and the temperature value corresponding to the corresponding pixel in the first temperature map exceeds a preset temperature threshold. If so, correct the temperature value corresponding to the corresponding pixel in the second temperature map.

[0022] The steps for performing magnetic susceptibility correction on the phase difference diagram include:

[0023] Determine whether the absolute value of the difference between the phase difference value corresponding to each pixel in the calibrated phase difference map and the phase difference value corresponding to the corresponding pixel in the reference phase map exceeds a preset phase difference threshold. If so, correct the phase difference corresponding to the corresponding pixel in the calibrated phase difference map.

[0024] Optionally, the ablation calculation method of the present invention further includes a step of correcting phase errors caused by motion, which removes the phase errors caused by motion by using linear least squares fitting of at least two sets of phase maps corresponding to different echo times at each pixel.

[0025] Optionally, the ablation calculation method of the present invention further includes the step of obtaining a weighted temperature map, which is obtained by weighting a temperature map obtained from at least one reference phase difference map and a temperature map obtained from at least one calibrated phase difference map. The echo time corresponding to the at least one calibrated phase difference map is not less than 20ms, 19ms, 18ms, 17ms, 16ms, 15ms, 14ms, 13ms, or 12ms.

[0026] The above optional steps for eliminating phase drift caused by the magnetic resonance system, correcting errors caused by magnetic susceptibility, correcting phase errors caused by motion, and obtaining the weighted temperature map, whether selected individually, partially, or in whole as part of the ablation calculation method of the present invention, all fall within the scope of the present invention.

[0027] Optionally, in the ablation calculation method of the present invention, the ablation of pixels is calculated using the following formula:

[0028]

[0029] Among them, E a The activation energy is represented by A, the frequency factor is R, the universal gas constant is T(τ), the temperature (°C) is a function of time τ, t is the current time, and pixels with an Ω value exceeding a set threshold (e.g., 1) are considered ablated. Other methods and parameters for temperature-based ablation calculations known to those skilled in the art can also be used as alternatives and are part of this invention.

[0030] In a second aspect of the invention, a storage medium is provided, on which program code is stored, which, when executed, implements the ablation calculation method of the invention.

[0031] In a third aspect of the invention, an ablation calculation system is provided, comprising an ablation calculation module capable of executing the ablation calculation method of the present invention. Further, the ablation calculation module may further comprise:

[0032] The information acquisition module is used to receive or acquire magnetic resonance information, which includes at least a phase map corresponding to the echo time obtained by scanning the target under test using a gradient echo sequence containing i different echo times, where i is a positive integer greater than or equal to 2.

[0033] The temperature difference calculation module is used to select at least two sets of phase diagrams corresponding to different echo times to obtain the corresponding temperature difference diagram.

[0034] A temperature map calculation module is used to obtain a temperature map based on the temperature difference map, which may be a weighted temperature map.

[0035] The ablation calculation module is used to calculate the ablation status of each pixel based on the temperature map.

[0036] A phase drift correction module, which performs the steps to eliminate phase drift caused by the magnetic resonance system.

[0037] The magnetic susceptibility error correction module is used to perform steps to correct errors caused by magnetic susceptibility.

[0038] Motion error correction module, which performs the steps to correct phase errors caused by motion.

[0039] In a fourth aspect of the invention, a laser interstitial hyperthermia device is provided, comprising: a host, a laser ablation device, and an optical fiber assembly. The host includes a processor loaded with program code, which, when executed, implements the ablation calculation method of the present invention.

[0040] The innovative aspects of the embodiments of the present invention include one or more of the following:

[0041] 1. Eliminated the problem of errors occurring when wrapping the phase difference diagram of echo sequences corresponding to some echo times;

[0042] 2. Eliminated the problem of abnormal temperature and inability to display temperature and melting status due to abnormal magnetic susceptibility in some areas caused by heating;

[0043] 3. It eliminates abnormal temperature errors caused by the movement of cerebrospinal fluid, heartbeat, etc., such as abnormal temperature in the ventricles.

[0044] 4. The method of the present invention uses a non-tracing algorithm or an iterative algorithm, which has a small amount of computation, saves calculation time, and can quickly obtain temperature and ablation results. Attached Figure Description

[0045] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of this application. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0046] Figure 1 The amplitude, phase, and temperature maps of magnetic resonance data obtained by existing technologies in in vitro and in vivo environments are presented.

[0047] Figure 2 A schematic flowchart of an ablation calculation method provided for one embodiment of this application;

[0048] Figure 3 The phase map, phase difference map, and temperature map are provided for one embodiment of this application;

[0049] Figure 4 A partial flowchart illustrating an ablation calculation method provided for another embodiment of this application;

[0050] Figure 5 A schematic diagram of an experimental apparatus provided for one embodiment of this application;

[0051] Figure 6 A partial magnification of the laser interstitial hyperthermia temperature map of an in vitro pork experiment provided as an embodiment of this application;

[0052] Figure 7 A tissue mimicry provided for one embodiment of this application ( Figure 7 (a) or extracted pork ( Figure 7 (b) Schematic diagram of temperature changes over time during the experiment;

[0053] Figure 8 A representative temperature graph of a dog (01) in an in vivo experiment provided for one embodiment of this application;

[0054] Figure 9 Laser interstitial hyperthermia ablation results for three representative dogs (dogs 01-03) provided as an embodiment of this application;

[0055] Figure 10 A schematic diagram of ablation estimates predicted by the ablation calculation method for the ablation T2w image provided for one embodiment of this application and for all other six different ablation laser doses (dog 04-09). Detailed Implementation

[0056] Magnetic resonance thermography (MRT) can guide various energy delivery-based therapeutic approaches, such as laser interstitial hyperthermia, focused ultrasound (FIT), and radiofrequency ablation, monitoring the temperature of the target tissue and the ablation effect. This invention uses MRT-guided laser interstitial hyperthermia as an example to illustrate the method of this invention. MRT-guided laser interstitial hyperthermia is a minimally invasive treatment that creates new options for treating tumors located in surgically challenging sites (anatomical or functional). This method rapidly coagulates tissue and induces tumor cell necrosis through protein denaturation by applying temperatures of 50–80°C or higher for tens of seconds. Compared to open surgery, laser interstitial hyperthermia can more precisely target tumors, reduce discomfort and infection risks, and shorten hospital stays. During the ablation process of laser interstitial hyperthermia, concurrent MRT plays an important role in more effectively ablating tumor cells and better protecting healthy surrounding cells and critical structures. Most laser interstitial hyperthermia ablation procedures rely on thermometry based on proton resonance frequency displacement.

[0057] However, as described in the background section, existing magnetic resonance temperature imaging methods are greatly affected by environmental factors, which can easily lead to a large discrepancy between the final temperature map and the actual temperature.

[0058] During laser-assisted interstitial hyperthermia (LAH) ablation, the inventors discovered that the main sources of error in the acquired temperature maps are phase errors caused by unwinding misalignment, magnetic susceptibility errors, and phase errors caused by motion. As the ablation laser dose changes, the magnetic susceptibility leads to a decrease in image amplitude and corresponding errors in image phase, thereby disrupting the reconstructed temperature map of the heating center and its surroundings. Errors in the reconstructed temperature map can lead to inaccurate estimation of the ablation area, potentially resulting in variations in treatment efficacy and thermal damage to critical tissues. Therefore, accurate temperature imaging is crucial for the effectiveness and safety of LHA, especially when applied to tight ablation areas in brain tissue.

[0059] Further research by the inventors revealed that the thermometry method based on proton resonance frequency shift is grounded in the fact that the resonance frequency of hydrogen protons varies with temperature within water molecules. For aqueous tissues, the change in the local magnetic field with temperature can be described as follows:

[0060]

[0061] Where α is the proton resonance frequency coefficient as a function of temperature, taken here as 0.008-0.015 ppm / ℃. The corresponding change in the resonance frequency of water protons affected by temperature can be expressed as:

[0062] Δf=αγB0·ΔT; (2)

[0063] Where ΔT represents temperature change, Δf represents resonant frequency change, γ represents gyromagnetic ratio, and B0 represents static magnetic field strength.

[0064] Changes in resonant frequencies due to temperature variations can be observed in the phase of complex magnetic resonance imaging. For a given gradient echo sequence interval TE, the relative temperature change ΔT can be calculated from the phase difference Δφ, and this equation can be expressed as:

[0065]

[0066] Gradient echo sequences are sequences used in thermometry based on proton resonance frequency displacement. According to formula (3), the longer the gradient echo sequence, the greater the phase difference that the same temperature change may result in, indicating that higher temperature sensitivity can be obtained.

[0067] refer to Figure 1 As the echo time of the gradient echo sequence increases, both phase contrast and phase wrapping increase, indicating that temperature sensitivity is higher and phase unwrapping procedures are more frequent in later echo times. Figure 1 In (a) (ex vivo, pig brain) and (b) (in vivo), the amplitude (top row) and phase (second row) of the first to fourth echoes obtained using gradient echo sequences with four different echo times, as used in embodiments of this application, are shown. Temperature maps (bottom row) are calculated based on each TE (echo time) using a conventional PRF algorithm. More phase wrapping occurs at longer echo times, as the image contrast increases accordingly. Note that intense laser heating can lead to signal loss due to changes in magnetic susceptibility and can also translate into phase and temperature errors in pixels around the heating center. In in vivo experiments, note that cerebrospinal fluid (CSF) motion can cause unsuitable high temperatures on the MRTI, which is more pronounced in the previous echo because shorter TEs are less tolerant of the introduced similar phase errors.

[0068] For example, inter-scan motion can be a major problem in temperature maps measured using proton resonance frequency displacement thermometry due to the movement of cerebrospinal fluid (CSF) in the brain. The size and phase signals of CSF are frequently altered on pulse gradient echo sequences by the normal dynamic movement of CSF, which can confound temperature estimations. CSF motion can also cause pixel shifts within and around the ventricles, leading to errors in phase contrast maps. Figure 1 As shown in b, the body temperature map shows a pseudo-high temperature in the third ventricle due to cerebrospinal fluid movement. The temperature error is more pronounced on pulse gradient echo sequences with shorter echo times because they are less tolerant of the phase shift intensity introduced by cerebrospinal fluid flow in (3).

[0069] In fact, the local magnetic field of water protons should also take into account the magnetic susceptibility x0, and formula (1) becomes:

[0070]

[0071] Where, σ χ0 This indicates the change in the local magnetic field caused by magnetic susceptibility.

[0072] Further research revealed that laser heating causes significant magnetization artifacts in GRE imaging around the laser tip. (Still referencing...) Figure 1 Heating centers with rapidly changing temperatures (such as...) Figure 1 (A) As indicated by the arrow, severe signal loss is observed over a long echo time. The spin phase shift within the voxel is caused by local magnetic field inhomogeneities resulting from variations in temperature and magnetic susceptibility.

[0073] Magnetization artifacts caused by laser heating, especially those in images corresponding to gradient echo sequences with longer echo times, are a significant cause of error. (Still referencing...) Figure 1 In both in vitro and in vivo experiments, phase errors around the heating center translate into pseudo-cryogenics in magnetic resonance thermography. Typically, during the imaging process, it is recommended to use gradient pulse sequences with the shortest possible echo times to minimize susceptibility artifacts. However, gradient pulse sequences with longer echo times can offer better temperature sensitivity and signal-to-noise ratio, presenting a current dilemma.

[0074] To balance temperature sensitivity, signal-to-noise ratio, and low error, embodiments of this application provide an ablation calculation method, which includes:

[0075] The target under test is scanned using a gradient echo sequence containing i different echo times to obtain a phase map corresponding to the echo times, where i is a positive integer greater than or equal to 2;

[0076] Select at least two sets of phase diagrams corresponding to different echo times to obtain the corresponding temperature difference diagrams;

[0077] A temperature map is obtained based on the temperature difference map;

[0078] The ablation status of each pixel is calculated based on the temperature map.

[0079] The ablation calculation method obtains i sets of phase maps based on a gradient echo sequence containing i different echo times, selects at least two sets of phase maps corresponding to different echo times to obtain corresponding phase difference maps, and obtains a temperature map based on the temperature difference map. The inventors discovered that the echo time of the gradient echo sequence is directly proportional to the magnitude of the magnetic susceptibility artifact. Therefore, the phase map obtained from the gradient echo sequence corresponding to the smaller echo time is least affected by the change in magnetic susceptibility due to heating, and its image data still maintains the correct phase. Therefore, the temperature map can be obtained based on the i sets of phase maps and the phase difference map obtained from the gradient echo sequence containing i different echo times, thereby reducing the error of the final obtained temperature map and improving its accuracy.

[0080] Furthermore, the ablation calculation method is neither a retrospective nor an iterative algorithm, has a small computational load, can provide almost real-time temperature maps, and has high reference value.

[0081] This application provides an ablation calculation method, such as... Figure 2 As shown, it includes:

[0082] S101: Use a gradient echo sequence containing i different echo times to scan the target under test, and obtain i sets of phase maps corresponding to the echo times, where i is a positive integer greater than or equal to 2;

[0083] In step S101, the minimum and maximum values ​​of the echo time in the gradient echo sequence can be determined according to actual needs. Generally, to minimize magnetic susceptibility artifacts caused by changes in magnetic susceptibility, the minimum value of the echo time in the gradient echo sequence can be the minimum value that the magnetic resonance thermal imaging device can obtain. The maximum value of the echo time in the gradient echo sequence generally does not exceed the upper limit of the range of echo time values ​​for imaging the target under test. For example, for head imaging, the range of echo time values ​​for the selectable gradient time sequence is 3 to 30 ms, and the specific values ​​of the echo time included in the gradient echo sequence must be within this range.

[0084] refer to Figure 3 For details, please refer to Figure 3 (a) The obtained phase map is as follows Figure 3 As shown in (a). Then as Figure 3 As shown in (b), the phase difference diagram is calculated using a complex subtraction procedure. The complex subtraction method avoids problematic phase wrapping.

[0085] The gradient echo sequence information containing i different echo times can be obtained by reading or receiving from a server or other storage device, or it can be obtained in real time according to the settings of the staff. This application does not limit the specific method of obtaining the gradient echo sequence, but depends on the actual situation.

[0086] S102: Select at least two sets of phase diagrams corresponding to different echo times to obtain the corresponding temperature difference diagram;

[0087] Optionally, between steps S102 and S103, static magnetic field strength drift correction of the phase diagram and phase difference diagram may also be included to eliminate errors caused by static magnetic field strength.

[0088] S103: Obtain a temperature map based on the temperature difference map.

[0089] S104: Calculate the ablation status of each pixel based on the temperature map. Ablation can be calculated using various methods based on temperature and time parameters, such as using the following formula:

[0090]

[0091] Among them, E a Ω represents the activation energy, A is the frequency factor, R is the universal gas constant, T(τ) is the function of temperature (°C) and time τ, t is the current time, and pixels with an Ω value exceeding a set threshold (e.g., 1) are considered to have been ablated.

[0092] The following describes the feasible execution methods of each step of the ablation calculation method provided in the embodiments of this application.

[0093] Based on the above embodiments, in one embodiment of this application, the specific steps for obtaining a temperature difference map include:

[0094] The phase difference map for any given moment is obtained by subtracting the phase map of a reference moment from the phase map of any of the different moments. The reference moment is any moment before energy (e.g., thermal energy, optical energy, radiofrequency ablation, cryoablation) is transferred to the target tissue, preferably a moment shortly before energy transfer, such as the moment when energy transfer is about to take place. The phase difference map corresponding to at least one echo time at that moment is selected as the reference phase difference map, and the phase difference maps corresponding to other echo times are calibrated to obtain the calibrated phase difference map. The echo time corresponding to the reference phase difference map is less than the echo time corresponding to the calibrated phase difference map.

[0095] The temperature difference map at that moment is calculated using the reference phase difference map and the calibrated phase difference map.

[0096] Optionally, the echo time of the reference phase difference diagram is less than or equal to 18ms. Preferably, the echo time of the reference phase difference diagram is no more than 17ms, 16ms, 15ms, 14ms, 13ms, 12ms, 11ms, 10ms, 9ms, 8ms or 7ms. More preferably, the echo time of the reference phase difference diagram is no more than 6ms, 5ms or 4ms.

[0097] Optionally, the phase map corresponding to the minimum echo time in the gradient echo sequence can be used as a reference phase map, and a reference phase difference map can be obtained based on the reference phase map to minimize the influence of magnetic susceptibility changes caused by heating on the phase map. As mentioned above, the minimum echo time in the gradient echo sequence can be the minimum value that the magnetic resonance thermal imaging device can obtain.

[0098] The step of selecting the phase difference map corresponding to at least one echo time at that moment as a reference phase difference map and calibrating the phase difference maps corresponding to other echo times includes the following steps:

[0099] Using the reference phase difference map and the phase difference map corresponding to the echo time to be calibrated, based on the proportional relationship between phase difference and echo time, an estimated value of the phase difference of the phase difference map corresponding to the echo time to be calibrated is calculated based on the phase difference of the echo time and the reference phase difference map. Then, using the estimated value, the phase difference to be calibrated is deconvolved according to the phase periodicity to obtain the calibrated phase difference.

[0100] Based on the above embodiments, in another embodiment of this application, the ablation calculation method further includes:

[0101] S105: A step to eliminate phase drift caused by the magnetic resonance system, wherein the step to eliminate phase drift caused by the magnetic resonance system is performed on a phase difference map or a temperature difference map.

[0102] The steps for eliminating phase drift caused by a magnetic resonance system on a phase difference map include:

[0103] Multiple thermal reference points (ROIs) are selected, and the average phase difference of the thermal reference points is subtracted from each phase difference map;

[0104] The steps to eliminate phase drift caused by a magnetic resonance system on a temperature difference map include:

[0105] Correction is performed by subtracting the average temperature difference of any of the aforementioned thermal reference points from the temperature difference graph.

[0106] Based on the above embodiments, in another embodiment of this application, the ablation calculation method further includes:

[0107] S106: Steps for magnetic susceptibility correction, which are performed on a phase difference diagram or temperature difference diagram;

[0108] The steps for performing magnetic susceptibility correction on a temperature difference graph include:

[0109] A first temperature map is obtained based on the reference phase difference map, and a corresponding second temperature map is obtained based on the calibrated phase difference map.

[0110] Optionally, a preset region is determined in the first temperature map and each of the second temperature maps;

[0111] Determine whether the absolute value of the difference between the temperature value corresponding to each pixel in the second temperature map and the temperature value corresponding to the corresponding pixel in the preset area of ​​the first temperature map exceeds a preset temperature threshold. If so, correct the temperature value corresponding to the corresponding pixel in the second temperature map. The correction can be done in various ways, such as replacing the temperature value of the second temperature map with the temperature value of the first temperature map in the pixel, replacing the temperature value of the pixel in the second temperature map with the temperature value of the adjacent pixel in the second temperature map, or fitting an approximate temperature based on the temperature values ​​of the adjacent pixels and the temperature value of the first temperature map to replace the temperature value of the second temperature map.

[0112] The steps for performing magnetic susceptibility correction on the phase difference diagram include:

[0113] A preset region is determined in the reference phase difference map and the calibrated phase difference map;

[0114] Determine whether the absolute value of the difference between the phase difference value corresponding to each pixel in the calibrated phase difference map and the phase difference value corresponding to the corresponding pixel in the preset area of ​​the reference phase map exceeds the preset phase difference threshold. If so, the phase difference in the calibrated phase difference map is corrected. The correction method is similar to that described above and will not be repeated.

[0115] Based on the above embodiments, in another embodiment of this application, the ablation calculation method further includes:

[0116] S107: Steps for correcting phase errors caused by motion on a phase difference diagram or temperature difference diagram;

[0117] The steps for correcting phase errors caused by motion on the phase difference map include:

[0118] Motion-induced phase errors are removed by using linear least-squares fitting of the reference phase difference map and the calibrated phase difference map at each pixel;

[0119] The steps for correcting phase errors caused by motion on a temperature difference map include:

[0120] A first temperature map is obtained based on the reference phase difference map, and a corresponding second temperature map is obtained based on the calibrated phase difference map.

[0121] Motion-induced phase errors are removed by using linear least-squares fitting of the first and second temperature maps at each pixel.

[0122] Still referencing Figure 4 (c), Figure 4 (c) shows the phase difference (first row) and relative temperature change (second row) as a function of time with and without motion error correction (left panel) and with (right panel). For shorter echo times, the phase error Δφ(x, y) bias It will introduce a large temperature deviation, but it can be correctly eliminated after linear least squares fitting.

[0123] As mentioned above, steps S105, S106, and S107 can be performed on either the phase difference layer or the temperature difference map. That is, the step of eliminating phase drift caused by the magnetic resonance system is performed on the phase difference map and / or the temperature difference map, the step of magnetic susceptibility correction is performed on the phase difference map and / or the temperature difference map, and the step of correcting phase errors caused by motion is performed on the phase difference map and / or the temperature difference map.

[0124] Still referencing Figure 3 , Figure 3 (c) shows the phase difference diagram after unwrapping and drift correction. Figure 3 (d) represents the temperature graph. Figure 3 (e) represents the image after magnetic susceptibility correction. Figure 3 (f) represents the image after motion error correction.

[0125] Based on the above embodiments, in another embodiment of this application, obtaining the temperature map based on the temperature difference map includes:

[0126] S1031: Calculate the temperature using the reference phase difference map and the calibrated phase difference map, and weight the calculated temperatures to obtain a temperature map of the target under test;

[0127] or

[0128] The average temperature difference is obtained by weighting the reference phase difference map and the calibrated phase difference map, and the temperature map of the target under test is calculated based on the average temperature difference.

[0129] In step S1031, the weighting can be various weighting methods, such as average weighting, or the temperature map can be a temperature map corresponding to a single echo time, that is, the weighting coefficient of the temperature map is 1, and the weighting coefficient of the temperature maps of other phases is 0.

[0130] The following steps may also be included after step S1031:

[0131] S108: Perform multiple interpolation processes on the temperature map of the target to be tested, and use the interpolated temperature map of the target to be tested to calculate the boundary of the ablation region.

[0132] The purpose of performing multiple interpolation processes on the temperature map of the target to be measured is to obtain a smoother ablation region boundary. The specific number of interpolation processes can be 2 or 3.

[0133] The following formula is used specifically in the calculation of the ablation region boundary:

[0134]

[0135] Among them, E a Ω represents the activation energy, A is the frequency factor, R is the universal gas constant, T(τ) is the temperature (°C) as a function of time τ, and t is the current time. Pixels with an Ω value exceeding a set threshold (e.g., 1) are considered to have been ablated.

[0136] The ablation calculation method provided in the embodiments of this application will be verified by specific experiments below.

[0137] The gel tissue mimic (gel phantom) was heated using a laser ablation system comprising a 10W, 980nm diode laser and a cooled laser applicator system. Phase images were acquired using a 3T MR scanner (Ingenia, Philips Healthcare, Best, Netherlands) with 16 receiver coils and a multi-echo temporal gradient echo sequence: flip angle = 30°, TE = 6 / 12 / 18 / 24ms, TR = 22ms, matrix = 176×176, FOV = 200x200mm. 2 Slice thickness = 5mm, 3s / image.

[0138] like Figure 5 As shown, two MR-compatible fiber optic temperature probes were also inserted into the tissue simulator, with the probe tips positioned near the ablation fiber to obtain the gel temperature at various points. Because the fiber optic probes were affected by the ablation fiber during heating, the thermometers only monitored the cooling phase. Figure 5 middle. Specifically, Figure 5 The image shows that an ablation fiber and two fiber optic temperature probes were inserted into the gel tissue simulator, and a reference tube filled with gel was fixed around it as an insulating reference.

[0139] In vitro experiments on pork and pig brain were conducted using the same scanning parameters as those used in tissue simulant experiments. Two experiments were performed for each tissue type (gel, pork, pig brain), one involving several laser cycles of heating and the other involving continuous heating and cooling. The root mean square error between the MR measurement temperature and the fiber optic measurement temperature was calculated as the measurement value for temperature accuracy.

[0140] The in vivo experiments on Doberman Pinschers have been approved by the Ethics Review Committee of Tsinghua University. Nine adult Doberman Pinschers received laser interstitial hyperthermia. The heating process was monitored using a multi-echo time gradient echo sequence with 32 receiver head coils on a 3T MR scanner (Ingenia, Philips Healthcare, Best, Netherlands).

[0141] Following the ablation procedure, post-ablation images were acquired to obtain detailed information about the actual extent of the ablation zone, including T1-gadolinium (T1+Gd) contrast images, fluid attenuation inversion recovery (FLAIR) images, diffusion-weighted MR (distorted, corrected by supplementation with FSL 5.0 or by EPSI method), and T2-weighted images.

[0142] Still referencing Figure 3 and Figure 4 , Figure 3 and Figure 4 This paper illustrates an example of temperature calculation using the ablation calculation method provided in the embodiments of this application. Figure 3 (a) During laser thermotherapy, a time frame was obtained by first acquiring a coil combination phase image via a multi-TE echo sequence; Figure 3 (b) Then a phase difference map is obtained, with white arrows indicating phase wrapping that occurs on the phase map around the heating center; Figure 3 (c) Displays the phase difference map after phase unwrapping and B0 drift correction. Figure 3 (d) According to the PRF offset method from Figure 3 (c) Calculated temperature plot. The white arrows highlight the error caused by magnetic susceptibility. Figure 3 (e) Temperature plot after magnetic susceptibility correction. The white arrows indicate errors caused by residual CSF motion. Figure 3 (f) Temperature map of motion correction.

[0143] Figure 4An exemplary method flow for representative pixels includes: Step 1, acquiring a phase difference map and a phase map obtained by unwrapping a reference phase difference map (TE1), as shown by the black arrows. Under rapid temperature changes, the phase difference map corresponding to some echo times becomes wrapped. Step 2, acquiring the unwrapped phase map. Step 3, static magnetic field strength (magnetic resonance system) drift correction, i.e., B0 drift correction, is to reduce system fluctuations. Step 4), phase error correction caused by magnetic susceptibility, using the shortest echo time (TE) to correct the temperature error caused by changes in magnetic susceptibility (black arrows) over longer echo times. Step 5, phase error correction caused by motion. The graphs in the first and second rows show the phase difference and corresponding temperature changes over time, respectively. For multiple echo times, the phase error caused by motion (black arrows) is almost the same, thus leading to a more pronounced temperature error over shorter TEs. The result of correcting motion errors shows smoother phase and temperature curves.

[0144] Tissue simulants and in vitro experimental results:

[0145] Figure 6 Representative temperature maps of ex vivo pork samples during laser interstitial hyperthermia are shown. Six representative images were selected from 300 frames (3s / frame) acquired during the thermal cycling (#50 represents frame 50, #146 represents frame 146, and so on). The first and second rows show the results using a conventional phase unwrapping method and a multi-echo time-based phase unwrapping method (multi-TE unwrapping) proposed in this embodiment, respectively. Using the existing phase unwrapping method, pixels on the temperature map are severely damaged due to magnetic susceptibility changes caused by laser heat, and cannot be recovered even after laser use. The technical principle is as follows: the existing phase unwrapping method is applied to the time dimension for phase jump detection. If the phase difference map of the current frame is incorrectly unwrapped, all subsequent frames will be affected. On the other hand, the ablation calculation method proposed in this invention is based on a multi-echo dimension, thus avoiding interference from previous frames. The third row shows a single echo time temperature map after phase unwrapping and magnetic susceptibility correction, with damaged pixels around the heating center correctly recovered. The last row shows the result of combining multi-echo time data using the ablation calculation method provided in this embodiment. The final magnetic resonance thermography showed a more uniform temperature at the hot spot.

[0146] Figure 7 The gel tissue mimicry (shown) Figure 7 (a) or extracted pork ( Figure 7 (b) The temperature changes over time during the experiment, measured by two temperature-measuring optical fibers (red lines) and calculated by the method provided in this application (dashed black lines). During multiple heating cycles ( Figure 7 (a)) or single heating ( Figure 7 In case (b), the temperature-time behavior calculated by the proton resonance frequency (PRF) closely matches the measured temperature-time behavior during the cooling phase with that measured by the fiber optic temperature probe (also known as the temperature-sensing fiber). Table 2 lists the root mean square error (RMSE) values ​​between the calculated MR values ​​and the temperature-sensing fiber measurements, representing the temperature accuracy of the proposed algorithm. Experiment 1 involved several laser cycles of heating, while Experiment 2 consisted of continuous heating and cooling phases. The results show that, in most cases, the RMSE error for the gel, pork, or pig brain tissue was less than 0.5 °C. Figure 7 In the image, the probe (left) represents the left-side temperature-sensing fiber optic cable, and the probe (right) represents the right-side temperature-sensing fiber optic cable.

[0147] Table 2. Comparison between the temperature measured by the optical fiber and the temperature calculated by MR using the method provided in the embodiments of this application.

[0148]

[0149] Abbreviation: RMSE, root mean square error; experiment, experiment L(R), left (right) fiber optic temperature probe.

[0150] Figure 8 Representative temperature maps of dog 01 in in vivo experiments are shown. It is important to note that the ablation area is located near the third and lateral ventricles. 100 frames (3 s / frame) acquired during laser ablation are superimposed on a post-ablation T2w MRI. From top to bottom, the temperature maps are calculated using existing algorithms based on single echo time (TE) data (TE = 6 ms and TE = 24 ms), and the temperature maps calculated using the algorithm proposed in this invention based on multiple (joint) TE echo sequences. The first row (TE = 6 ms) shows pseudo-high temperatures within the third and lateral ventricles, indicating that the short TE calculated temperature is severely affected by CSF flow artifacts. CSF-induced artifacts within the third ventricle (indicated by white arrows) are still present in the second row (TE = 24 ms), but are well suppressed by the proposed joint TE echo sequence algorithm. The second line shows that a longer TE (TE = 24 ms) provides smoother boundaries and better temperature SNR compared to a shorter TE (TE = 6 ms). However, as mentioned above, pixels around the heating center are damaged due to variations in magnetic susceptibility. On the other hand, our proposed method integrates information from multiple echoes, thus the obtained temperature map simultaneously eliminates errors caused by CSF and magnetic susceptibility, revealing a more uniform and symmetrical heating region.

[0151] Figure 9The results of laser interstitial thermotherapy ablation in three representative dogs (dogs 01-03) are shown. The first column is the T2w image (post-T2) after thermal ablation and gadolinium application. The location of the ablation fiber is clearly shown on the T2w MRI. The second column shows the final estimated ablation damage at a given laser dose using the method proposed in this invention, which is superimposed in red (darker color on the grayscale image) on the post-ablation T2w image. This demonstrates good agreement between the estimated ablation calculated by the method proposed in this invention and the post-ablation MRI. The ablation area estimated by the method of this invention is shown in the upper left corner of the ablation image. The third column is a representative temperature map obtained when laser heating is most intense. The last three columns are FLAIR, DWI, and T1w images after ablation with gadolinium application, respectively. They all show a sharp transition zone between dead and living tissue.

[0152] Figure 10 The image shows post-ablation T2w images and the final algorithm-predicted damage estimates for all six other cases (dog 04-09) with different ablation laser doses. The calculated ablation values ​​(second row) match well with the post-ablation assessment values ​​(first row). The ablation area estimated by the algorithm is shown in the upper left corner of the T2w image. Depending on the duration of laser heating, the laser ablation area ranges from less than 30 square millimeters to nearly 90 square millimeters.

[0153] The experimental results above demonstrate that the ablation calculation method provided in this application can be used to correct errors caused by the magnetic susceptibility in the proton resonance frequency temperature map due to the heating laser itself. We first propose applying multi-echo time-gradient echo pulse sequences to magnetic resonance thermography using a proton resonance frequency shifting method. Compared to single-echo sequences, multi-gradient echo sequences provide more information without additional scan time and offer novel methods for phase unrolling and artifact removal.

[0154] Shorter echo times can tolerate magnetic susceptibility artifacts but are sensitive to noise, while longer echo times offer better temperature sensitivity and signal-to-noise ratio but are significantly affected by magnetic susceptibility artifacts. The ablation calculation method proposed in this invention integrates the advantages of different echo times to obtain better temperature map measurement results. Furthermore, this ablation calculation method can significantly improve the robustness and signal-to-noise ratio of magnetic resonance thermography, thereby avoiding damage to healthy tissues due to miscalculation of low temperatures.

[0155] The method of this invention also exhibits excellent suppression of cerebrospinal fluid flow errors and can provide accurate temperature measurements within or around the ventricles. Compensating for errors caused by cerebrospinal fluid movement is clinically important for laser interstitial hyperthermia in the treatment of periventricular brain lesions. Furthermore, the proposed algorithm is online-compatible and does not require iterative computation, making it well-suited for magnetic resonance thermography, where very close to real-time temperature maps are required.

[0156] The magnetic resonance thermal imaging system provided in the embodiments of this application is described below. The magnetic resonance thermal imaging system described below can be referred to in conjunction with the ablation calculation method described above.

[0157] Accordingly, this application also provides an ablation calculation system, which includes an ablation calculation module capable of executing the ablation calculation method of the present invention. Further, the ablation calculation module may further include:

[0158] The information acquisition module is used to receive or acquire magnetic resonance information, which includes at least a phase map corresponding to the echo time obtained by scanning the target under test using a gradient echo sequence containing i different echo times, where i is a positive integer greater than or equal to 2.

[0159] The temperature difference calculation module is used to select at least two sets of phase diagrams corresponding to different echo times to obtain the corresponding temperature difference diagram.

[0160] A temperature map calculation module is used to obtain a temperature map based on the temperature difference map, which may be a weighted temperature map.

[0161] The ablation calculation module is used to calculate the ablation status of each pixel based on the temperature map.

[0162] A phase drift correction module, which performs the steps to eliminate phase drift caused by the magnetic resonance system.

[0163] The magnetic susceptibility error correction module is used to perform steps to correct errors caused by magnetic susceptibility.

[0164] Motion error correction module, which performs the steps to correct phase errors caused by motion.

[0165] In a fourth aspect of the invention, a laser interstitial hyperthermia device is provided, comprising: a host, a laser ablation device, and an optical fiber assembly. The host includes a processor loaded with program code, which, when executed, implements the ablation calculation method of the present invention.

[0166] This application also provides another magnetic resonance thermal imaging system, including:

[0167] The data transmission module is configured to receive magnetic resonance imaging sequences and determine image integrity.

[0168] The temperature calculation module is configured to select a sequence, calculate the phase difference, calibrate the phase difference, and calculate the temperature.

[0169] The temperature display module is configured to display the temperature in a pseudo-color graph or isotherm mode.

[0170] The ablation calculation module is configured to calculate and display the ablation results;

[0171] The system takes no more than 1 second to perform a complete calculation.

[0172] In some embodiments of this application, the time for the system to perform a complete calculation is preferably no more than 0.5s, and most preferably no more than 0.1s.

[0173] Accordingly, embodiments of this application also provide a magnetic resonance thermal imaging system, including: a memory and a processor;

[0174] The memory is used to store program code, and the processor is used to call the program code, which is used to execute the ablation calculation method described in any of the above embodiments.

[0175] Accordingly, this application also provides a storage medium storing program code, which, when executed, implements the ablation calculation method described in any of the above embodiments.

[0176] In summary, the embodiments of this application provide an ablation calculation method and related apparatus. The ablation calculation method is not a traceability algorithm or an iterative algorithm, has a small computational load, can provide almost real-time temperature maps, and has high reference value.

[0177] The features described in the various embodiments of this specification can be substituted for or combined with each other. Each embodiment focuses on the differences from other embodiments, and the same or similar parts between the various embodiments can be referred to each other.

[0178] The above description of the disclosed embodiments enables those skilled in the art to make or use this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. An ablation calculation method, characterized in that, include: The target under test is scanned using a gradient echo sequence containing i different echo times to obtain a phase map corresponding to the echo times, where i is a positive integer greater than or equal to 2; One echo time from the i different echo times of the gradient echo sequence is selected as the reference echo time. The value of the reference echo time is less than or equal to 18ms. The phase diagram corresponding to the reference echo time at any time of the gradient echo sequence is subtracted from the phase diagram corresponding to the reference echo time at the reference time to obtain a reference phase difference diagram. The phase difference diagrams corresponding to other echo times are calibrated based on the reference phase difference diagram to obtain a calibrated phase difference diagram. The temperature difference diagram is calculated using the reference phase difference diagram and the calibrated phase difference diagram. A temperature map is obtained based on the temperature difference map; The ablation process is calculated based on the temperature graph.

2. The ablation calculation method according to claim 1, characterized in that, The echo time corresponding to the reference phase difference map is less than the echo time corresponding to the calibrated phase difference map.

3. The ablation calculation method according to claim 2, characterized in that, The calibration is performed as follows: Based on the proportional relationship between the phase difference diagram and the echo time, the estimated value of the phase difference diagram to be calibrated is calculated based on the echo time and the reference phase difference diagram. Using the estimated value, the phase difference map to be calibrated is unwrapped according to the phase periodicity to obtain the calibrated phase difference map.

4. The ablation calculation method according to any one of claims 1 to 3, characterized in that, It also includes steps to eliminate phase drift caused by the magnetic resonance system.

5. The ablation calculation method according to claim 4, characterized in that, In the step of eliminating phase drift caused by the magnetic resonance system, several regions with stable physical temperature and uniform tissue are selected as thermal reference points. Phase drift correction is performed by subtracting the average phase difference of the thermal reference points from each phase difference image or subtracting the average temperature change of the thermal reference points from the temperature change image.

6. The ablation calculation method according to any one of claims 1 to 3, characterized in that, It also includes a step to correct errors caused by magnetic susceptibility, which is performed on a phase difference diagram or temperature difference diagram. This step includes: The steps for performing magnetic susceptibility correction on a temperature difference graph include: A first temperature map is obtained based on the reference phase difference map, and a corresponding second temperature map is obtained based on the calibrated phase difference map. Determine whether the absolute value of the difference between the temperature value corresponding to each pixel in the second temperature map and the temperature value corresponding to the corresponding pixel in the first temperature map exceeds a preset temperature threshold. If so, correct the temperature value corresponding to the corresponding pixel in the second temperature map. The steps for performing magnetic susceptibility correction on the phase difference diagram include: Determine whether the absolute value of the difference between the phase difference value corresponding to each pixel in the calibrated phase difference map and the phase difference value corresponding to the corresponding pixel in the reference phase map exceeds a preset phase difference threshold. If so, correct the phase difference corresponding to the corresponding pixel in the calibrated phase difference map.

7. The ablation calculation method according to any one of claims 1 to 3, characterized in that, It also includes a step to correct phase errors caused by motion, which removes the phase errors caused by motion by using linear least-squares fitting of at least two sets of phase maps corresponding to different echo times at each pixel.

8. The ablation calculation method according to any one of claims 1 to 3, characterized in that, The temperature map is a weighted temperature map obtained by weighting a temperature map obtained from at least one reference phase difference map and a temperature map obtained from at least one calibrated phase difference map.

9. The ablation calculation method according to claim 1, characterized in that, The ablation is calculated using the following formula: ； in, Indicates activation energy. It is the frequency factor, and R is the universal gas constant. It is temperature and time The function, where t is the current time. Pixels with values ​​exceeding a set threshold are considered to have been ablated.

10. A storage medium, characterized in that, The storage medium stores program code, which, when executed, implements the ablation calculation method according to any one of claims 1 to 9.

11. An ablation calculation system, characterized in that, It includes an ablation calculation module, which is capable of performing the ablation calculation method according to any one of claims 1 to 9.

12. A laser interstitial hyperthermia device, characterized in that, include: The device comprises a host computer, a laser ablation device, and an optical fiber assembly. The host computer includes a processor loaded with program code, which, when executed, implements the ablation calculation method according to any one of claims 1-9.

Images