CN114764133B - An ablation calculation method and ablation calculation system - Google Patents

Abstract

The present application discloses an ablation calculation method and an ablation calculation system. The ablation calculation method includes: using a gradient echo sequence containing i different echo times to scan the target to be measured, and obtaining a phase diagram corresponding to the echo time , i is a positive integer greater than or equal to 2; select at least two groups of phase diagrams corresponding to different echo times to obtain corresponding temperature difference diagrams; obtain a temperature diagram according to the temperature difference diagram; calculate the ablation situation according to the temperature diagram; ablation The calculation system includes an ablation calculation module capable of executing the ablation calculation method of the present invention.

Description

An ablation calculation method and ablation calculation system

technical field

The present application relates to the field of medical devices, and more specifically, to an ablation calculation method and an ablation calculation system.

Background technique

Magnetic Resonance Temperature Imaging (MRI) can realize non-invasive, real-time, in-vivo monitoring of temperature distribution and changes inside the object under test. It has important uses during ablation monitoring of therapy.

One of the current magnetic resonance temperature imaging methods is the temperature measurement method based on the displacement of the proton resonance frequency (Proton Resonance Frequency, PRF). The distribution of magnetic susceptibility is uneven, and objective environmental factors such as tissue movement caused by breathing/blood flow pulsation are greatly affected. Common errors caused include phase unwrapping errors, errors caused by rapid changes in magnetic susceptibility, and errors caused by movement. This causes a big difference between the finally obtained temperature map and the actual temperature, which makes the temperature map lose its reference value.

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.

Contents of the invention

In order to solve the above technical problems, the present application provides an ablation calculation method and a related ablation calculation system.

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

Scanning the target to be measured using a gradient echo sequence containing i different echo times to obtain a phase diagram corresponding to the echo time, where i is a positive integer greater than or equal to 2;

Selecting at least two groups of phase diagrams corresponding to different echo times to obtain corresponding temperature difference diagrams;

Obtaining a temperature map according to the temperature difference map;

The ablation (for each pixel) is calculated from the temperature map.

Further, in this method, the temperature difference map is obtained as follows: use the phase map at any time to subtract the phase map at the reference time to obtain the phase difference map at this moment, and select the phase difference map corresponding to at least one echo time as Refer to the phase difference diagram, based on the reference phase difference diagram, calibrate the phase difference diagram to be calibrated corresponding to other echo times to obtain a calibrated phase difference diagram, and the echo time corresponding to the reference phase difference diagram is smaller than the calibrated phase difference diagram The corresponding echo time; the temperature difference map at that moment is calculated using the reference phase difference map and the calibrated phase difference map.

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

Furthermore, in this method, the above-mentioned calibration is performed as follows:

According to the proportional relationship between the phase difference map and the echo time, the estimated value of the phase difference map to be calibrated is calculated based on the echo time and the reference phase difference map;

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.

Optionally, the ablation calculation method of the present invention also includes the step of eliminating the phase drift caused by the magnetic resonance system (for example, B0 drift error); further, in the step of eliminating the phase drift caused by the magnetic resonance system, select several physical temperature stable The area with no change and uniform tissue is used as a thermal reference point, and the phase drift is performed by subtracting the average phase difference of the thermal reference point from each phase difference image or subtracting the average temperature change of the thermal reference point from the temperature change map Correction.

Optionally, the ablation calculation method of the present invention further includes the step of correcting the error caused by the magnetic susceptibility, and the magnetic susceptibility correction step is performed on the phase difference map or the temperature difference map, and the step includes:

The steps to perform a magnetic susceptibility correction on a temperature difference map include:

Obtaining a first temperature map according to the reference phase difference map, obtaining a corresponding second temperature map according to the calibrated phase difference map,

judging 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; The temperature value corresponding to the corresponding pixel in the second temperature map is corrected;

The steps for susceptibility correction on a phase difference map include:

Judging 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, and if so, then The phase difference corresponding to the corresponding pixel in the calibration phase difference map is corrected.

Optionally, the ablation calculation method of the present invention further includes a step of correcting the phase error caused by motion, which step uses at least two sets of phase maps corresponding to different echo times at each pixel to fit the motion caused phase error removal.

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 at least one temperature map obtained from a calibrated phase difference map . Wherein, the echo time corresponding to at least one calibrated phase difference map is not less than 20ms, 19ms, 18ms, 17ms, 16ms, 15ms, 14ms, 13ms or 12ms.

The above optional step of eliminating the phase drift caused by the magnetic resonance system, the step of correcting the error caused by the magnetic susceptibility, the step of correcting the phase error caused by the motion, and the step of obtaining the weighted temperature map are individually, partially or completely selected as the present invention. The solutions of a part of the ablation calculation method all belong to the scope of the content of the present invention.

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

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 the value of Ω exceeds the set threshold (such as 1) Pixels are considered ablated. Other methods and parameters for temperature-based ablation calculations known to those skilled in the art may also be used as alternatives and form part of the present invention.

In the second aspect of the present invention, a storage medium is provided, on which a program code is stored, and when the program code is executed, the ablation calculation method of the present invention is realized.

In a third aspect of the present invention, an ablation calculation system is provided, which includes an ablation calculation module, and the ablation calculation module is capable of executing the ablation calculation method of the present invention. Further, the ablation calculation module may further include:

An information collection module, configured to receive or collect magnetic resonance information, the magnetic resonance information at least includes a phase map corresponding to the echo time obtained by scanning the target to be measured using a gradient echo sequence containing i different echo times, The i is a positive integer greater than or equal to 2;

A temperature difference calculation module, which is used to select at least two groups of phase diagrams corresponding to different echo times to obtain corresponding temperature difference diagrams;

A temperature map calculation module, which is used to obtain a temperature map according to the temperature difference map, and the temperature map may be a weighted temperature map;

An ablation calculation module, configured to calculate the ablation condition of each pixel according to the temperature map.

a phase drift correction module for performing the step of eliminating phase drift caused by the magnetic resonance system,

a susceptibility error correction module for performing the step of correcting susceptibility-induced errors,

A motion error correction module is configured to perform a step of correcting a phase error caused by motion.

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

The innovations of the embodiments of the present invention include one or more of the following:

1. Eliminated the problem that errors occurred when wrapping the phase difference diagram of the echo sequence corresponding to a part of the echo time;

2. Eliminate the problem that the abnormal magnetic susceptibility in some areas caused by the temperature rise causes abnormal temperature, and the temperature and ablation cannot be displayed;

3. Eliminate abnormal temperature errors caused by cerebrospinal fluid, heart beat and other movements, such as abnormal temperature in the ventricle.

4. The method of the present invention uses a non-retrospective algorithm or an iterative algorithm, which has a small amount of calculation, saves calculation time, and can quickly obtain temperature and ablation results.

Description of drawings

In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the following will briefly introduce the drawings that need to be used in the description of the embodiments or the prior art. Obviously, the accompanying drawings in the following description are only It is an embodiment of the present application, and those skilled in the art can also obtain other drawings according to the provided drawings without creative work.

Fig. 1 is the amplitude, phase diagram and temperature diagram of the magnetic resonance data obtained in the ex vivo and in vivo environment in the prior art;

FIG. 2 is a schematic flowchart of an ablation calculation method provided by an embodiment of the present application;

Fig. 3 is the obtained phase diagram, phase difference diagram and temperature diagram provided by one embodiment of the present application;

Fig. 4 is a schematic flowchart of a part of an ablation calculation method provided by another embodiment of the present application;

Fig. 5 is the schematic diagram of the experimental device provided by an embodiment of the present application;

Fig. 6 is the partial enlargement of the laser interstitial hyperthermia temperature map of the isolated pork experiment provided by one embodiment of the present application;

Figure 7 is a schematic diagram of the change of temperature over time during the experiment of gel tissue mimic (Figure 7(a)) or isolated pork (Figure 7(b)) provided by an embodiment of the present application;

Fig. 8 is a representative temperature diagram of dog 01 in an in vivo experiment provided by an embodiment of the present application;

Fig. 9 is the result of laser interstitial hyperthermia ablation of 3 representative dogs (dog 01-03) provided by one embodiment of the present application;

Fig. 10 is a schematic diagram of the ablation estimation predicted by the ablation calculation method of the T2w image after ablation and all other six different ablation laser doses (Dog 04-09) provided by an embodiment of the present application.

Detailed ways

Magnetic resonance temperature imaging can guide a variety of energy delivery treatments, such as laser interstitial hyperthermia, focused ultrasound therapy, radiofrequency ablation, etc., to monitor the target tissue temperature and ablation effect. The present invention uses magnetic resonance temperature imaging-guided laser interstitial hyperthermia as an example to illustrate the method of the present invention. Magnetic resonance temperature imaging-guided laser interstitial hyperthermia is a minimally invasive treatment method for the treatment of surgically challenging Tumors at different sites (anatomical or functional) create new options. This method can rapidly coagulate tissue and induce tumor cell necrosis through protein denaturation by applying a temperature of 50-80° C. or higher for tens of seconds. Compared with open surgery, laser interstitial hyperthermia targets tumors more precisely, reduces discomfort and infection risk, and shortens patient hospital stays. During the ablation process of laser interstitial hyperthermia, concurrent magnetic resonance thermography plays an important role for more effective ablation of tumor cells and better preservation of healthy surrounding cells and critical structures. Most laser interstitial hyperthermia ablation procedures depend on thermometry based on proton resonance frequency shifts.

However, as described in the background art, the existing magnetic resonance temperature imaging method is greatly affected by factors such as the environment, which may easily cause a problem that the finally obtained temperature map differs greatly from the actual temperature.

During the ablation process of laser interstitial hyperthermia, the inventor found through research that the main sources of errors in the obtained temperature map are the phase error caused by unwinding misalignment, the magnetic susceptibility error and the phase error caused by motion. As the ablative laser dose is changed, magnetic susceptibility leads to a reduction in image amplitude and a corresponding error in image phase, corrupting the reconstructed temperature map in and around the heating center. Errors in reconstructing the temperature map can lead to misestimation of the ablation area, which can lead to variations in treatment efficacy as well as thermal damage to critical tissues. Therefore, accurate temperature imaging is crucial for the effectiveness and safety of laser interstitial hyperthermia, especially when laser interstitial hyperthermia is applied to tight ablated regions in brain tissue.

The inventor found through further research that the thermometry method based on the shift of the proton resonance frequency is based on the fact that the resonance frequency of the hydrogen proton changes with the temperature in the water molecule. For hydrated tissues, the variation of the local magnetic field with temperature can be described as:

Wherein, α is the proton resonance frequency coefficient that varies with temperature, which is taken as 0.008-0.015ppm/°C here. The corresponding resonance frequency change of water protons affected by temperature can be expressed as:

Δf=αγB ₀ ·ΔT; (2)

Among them, ΔT represents the temperature change, Δf represents the resonance frequency change, γ represents the gyromagnetic ratio, and _B0 represents the static magnetic field strength.

Changes in resonance frequency due to temperature changes can be observed in the phase of complex magnetic resonance imaging. For a given interval time TE of the gradient echo sequence, the relative temperature change ΔT can be calculated according to the phase difference Δφ, and the equation can be expressed as:

A gradient echo sequence is a sequence used in thermometry based on the frequency shift of the proton resonance. According to formula (3), the longer the gradient echo sequence, the larger the phase difference may be caused by the same temperature change, indicating that higher temperature sensitivity can be obtained.

Referring to Fig. 1, as the echo time of the gradient echo sequence increases, both the phase contrast and the phase wrapping increase, indicating a higher temperature sensitivity and more phase unwrapping procedures at later echo times. In Fig. 1, in (a) (in vitro, pig brain) and (b) (in vivo), the first to fourth rounds obtained by the gradient echo sequence containing 4 different echo times used in the examples of this application Wave amplitude (top row) and phase (second row). The temperature map (lower row) is calculated for each TE (echo time) setting using the conventional PRF algorithm. More phase wrapping occurs at longer echo times because the image contrast increases accordingly. Note that intense laser heating causes signal loss due to changes in magnetic susceptibility, and also translates into phase and temperature errors in pixels around the heated center. In vivo experiments, note that Cerebrospinal Fluid (CSF) motion can cause unsuitably high temperature on MRI, which is more evident in the previous echo due to the similar phase error introduced by the shorter TE pair low tolerance.

For example, scan-to-scan motion can be a big problem in temperature maps measured by thermometry based on proton resonance frequency shifts due to the motion of cerebrospinal fluid in the brain. The magnitude and phase of the CSF signal are often altered on pulsed gradient-echo sequences by the normal dynamic movement of CSF, which may confound temperature estimates. CSF movement can also cause pixel shifts within and around the ventricles, resulting in errors in the phase contrast map. As shown in Figure 1b, the in vivo temperature map shows a pseudo-hyperthermia in the third ventricle due to CSF movement. Temperature errors are more pronounced on pulsed gradient-echo sequences with shorter echo times because they are less tolerant to the magnitude of the phase shift introduced by CSF flow in (3).

In fact, the local magnetic field of water protons should also consider the magnetic susceptibility x ₀ , the formula (1) becomes:

in, Indicates the local magnetic field change due to magnetic susceptibility.

Further investigation found that laser heating caused significant magnetization artifacts in GRE imaging around the laser tip. Still referring to FIG. 1 , heating centers with sharp temperature changes (as indicated by the arrows in FIG. 1( a )) show severe signal loss on the order of longer echo times. Intravoxel spin phase shifts are caused by local magnetic field inhomogeneities caused by temperature and magnetic susceptibility changes.

Magnetization artifacts caused by laser heating, especially in images corresponding to gradient echo sequences with longer echo times, are an important cause of errors. Still referring to Figure 1, in ex vivo or in vivo experiments, phase errors around the heating center translate into pseudo-cold temperatures on magnetic resonance thermography. Typically, during imaging, 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 provide better temperature sensitivity and signal-to-noise ratio, which is the current dilemma.

In order to take into account temperature sensitivity, signal-to-noise ratio, and low error, an embodiment of the present application provides an ablation calculation method, which includes:

Scanning the target to be measured using a gradient echo sequence containing i different echo times to obtain a phase diagram corresponding to the echo time, where i is a positive integer greater than or equal to 2;

Selecting at least two groups of phase diagrams corresponding to different echo times to obtain corresponding temperature difference diagrams;

Obtaining a temperature map according to the temperature difference map;

The ablation condition of each pixel is calculated according to the temperature map.

The ablation calculation method obtains i groups of phase images based on gradient echo sequences containing i different echo times, selects at least two groups of phase images corresponding to different echo times to obtain corresponding phase difference images, and according to the temperature difference images Get a temperature map. The inventor found through research that the echo time of the gradient echo sequence is proportional to the size of the magnetic susceptibility artifact, so the phase map obtained by the gradient echo sequence corresponding to the smaller echo time is magnetized due to heating The impact of rate changes is minimal, and its image data still maintains the correct phase. Therefore, the temperature map can be obtained based on the i group of phase maps and phase difference maps obtained from gradient echo sequences containing i different echo times, so as to reduce the final acquisition time. The error of the temperature map is used to improve the accuracy of the obtained temperature map.

Further, the ablation calculation method is not a retrospective algorithm or an iterative algorithm, and has a small amount of calculation, and can provide an almost real-time temperature map, which has high reference significance.

The embodiment of the present application provides an ablation calculation method, as shown in Figure 2, including:

S101: Scan the target to be measured using a gradient echo sequence containing i different echo times to obtain i groups of phase images corresponding to the echo times, where i is a positive integer greater than or equal to 2;

In step S101, the minimum value and the maximum value of the echo time in the gradient echo sequence can be determined according to actual needs. Generally, in order to minimize the magnetic susceptibility artifacts caused by the magnetic susceptibility change, the The minimum value of the echo time in the gradient echo sequence can take the minimum value that can be obtained by the magnetic resonance temperature imaging equipment, and the maximum value of the echo time in the gradient echo sequence generally does not exceed the echo of the imaging target to be measured. The upper limit of the time range. For example, for head imaging, the value range of the echo time of the optional gradient time series is 3-30 ms, and the specific values of the echo time included in the gradient echo sequence must be in this within the value range.

Referring to FIG. 3 , specifically referring to FIG. 3( a ), the obtained phase diagram is shown in FIG. 3( a ). Then, as shown in Fig. 3(b), the phase difference map is calculated by the complex subtraction procedure. Complex subtraction avoids problematic phase wrapping.

The acquisition of gradient echo sequence information containing i different echo times can be read or received from a server or other storage device, or can be obtained in real time according to the settings of the staff. The specific method is not limited and depends on the actual situation.

S102: Select at least two groups of phase diagrams corresponding to different echo times to obtain corresponding temperature difference diagrams;

Optionally, between step S102 and step S103, correction of static magnetic field strength drift of the phase map and phase difference map may also be included, so as to eliminate errors caused by the static magnetic field strength.

S103: Obtain a temperature map according to the temperature difference map.

S104: Calculate the ablation situation of each pixel according to the temperature map, the ablation can use various methods for calculating ablation based on temperature and time parameters, for example, use the following formula to calculate:

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 the value of Ω exceeds the set threshold (such as 1) Pixels are considered ablated.

A feasible implementation manner of each step of the ablation calculation method provided in the embodiment of the present application is described below.

On the basis of the above embodiments, in one embodiment of the present application, the specific steps of obtaining the temperature difference map include:

Use the phase diagram at any one of the different times to subtract the phase diagram at a reference time to obtain a phase difference diagram at that time, and the reference time is before energy is transmitted to the target tissue (such as thermal energy, light energy, radiofrequency ablation, cryoablation) At any time, preferably shortly before the energy transmission, such as the moment when the energy transmission is about to be carried out; at this moment, at least one phase difference diagram corresponding to the echo time is selected as the reference phase difference diagram, and the phase difference diagram corresponding to other echo times The phase difference diagram is calibrated to obtain a calibrated phase difference diagram, and the echo time corresponding to the reference phase difference diagram is smaller than the echo time corresponding to the calibrated phase difference diagram;

A temperature difference map at that moment is calculated using the reference phase difference map and the calibrated phase difference map.

Optionally, the value of the echo time of the reference phase difference diagram is less than or equal to 18ms, preferably, the value of the echo time of the reference phase difference diagram does not exceed 17ms, 16ms, 15ms, 14ms, 13ms, 12ms, 11ms, 10ms, 9ms, 8ms or 7ms, more preferably, the value of the echo time of the reference phase difference map is not more than 6ms, 5ms or 4ms.

Optionally, the phase map corresponding to the minimum echo time in the gradient echo sequence can be used as the reference phase map, and the reference phase difference map is obtained based on the reference phase map, so as to minimize The effect of magnetic susceptibility changes. As mentioned above, the minimum echo time in the gradient echo sequence may be the minimum value that can be obtained by the magnetic resonance temperature imaging equipment.

Selecting the phase difference diagram corresponding to at least one echo time at this moment as a reference phase difference diagram, and calibrating the phase difference diagram corresponding to other echo times includes the following steps:

Using the reference phase difference diagram and the phase difference diagram corresponding to the echo time to be calibrated, according to the proportional relationship between the phase difference and the echo time, based on the phase difference between the echo time and the reference phase difference diagram, the to-be-calibrated phase difference diagram is calculated to obtain the an estimated value of the phase difference of the phase difference diagram corresponding to the calibrated echo time; and then use the estimated value to deconvolute the phase difference to be calibrated according to the phase periodicity to obtain a calibrated phase difference.

On the basis of the above embodiments, in another embodiment of the present application, the ablation calculation method further includes:

S105: A step of eliminating the phase drift caused by the magnetic resonance system, the step of eliminating the phase drift caused by the magnetic resonance system is performed on a phase difference diagram or a temperature difference diagram.

The steps performed on the phase difference map to eliminate the phase drift caused by the magnetic resonance system include:

Select a plurality of thermal reference points (Region of Interest, ROI), by subtracting the average phase difference of the thermal reference points from each phase difference map;

Steps to remove phase drift induced by the MR system on the temperature difference map include:

Correction is performed by subtracting the average temperature difference at any one of said thermal reference points from the temperature difference map.

On the basis of the above embodiments, in another embodiment of the present application, the ablation calculation method further includes:

S106: the step of correcting the magnetic susceptibility, the step of correcting the magnetic susceptibility is performed on the phase difference diagram or the temperature difference diagram;

The steps to perform a magnetic susceptibility correction on a temperature difference map include:

obtaining a first temperature map according to the reference phase difference map, and obtaining a corresponding second temperature map according to the calibrated phase difference map;

Optionally, determining preset regions in the first temperature map and each of the second temperature maps;

Judging 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 in the first temperature map exceeds a preset temperature threshold, if so , then the temperature value corresponding to the corresponding pixel in the second temperature map is corrected; there are many methods for correction, for example, the temperature value of the first temperature map in the pixel can be used to replace the temperature value of the second temperature map, or use Replace 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 fit an approximate temperature based on the temperature value of the adjacent pixel and the temperature value of the first temperature map to replace the temperature value of the second temperature map temperature value;

The steps for susceptibility correction on a phase difference map include:

determining a preset region in the reference phase difference map and the calibrated phase difference map;

judging 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 in the reference phase map exceeds a preset phase difference threshold, If yes, correct the phase difference in the calibrated phase difference map, the correction method is similar to the above and will not be repeated.

On the basis of the above embodiments, in another embodiment of the present application, the ablation calculation method further includes:

S107: a step of correcting the phase error caused by motion on the phase difference map or the temperature difference map;

The steps to correct motion-induced phase errors on the phase difference map include:

removing motion-induced phase errors by a linear least squares fit at each pixel using the reference phase difference map and the calibrated phase difference map;

The steps to correct motion-induced phase errors on the temperature difference map include:

obtaining a first temperature map according to the reference phase difference map, and obtaining a corresponding second temperature map according to the calibrated phase difference map;

Motion-induced phase errors are removed by a linear least squares fit at each pixel using the first temperature map and the second temperature map.

Still referring to Figure 4(c), Figure 4(c) shows the phase difference (first row) and relative temperature change ( second line). For shorter echo times, the phase error Δφ(x,y) _bias introduces larger temperature deviations, but is correctly removed after linear least squares fitting.

As mentioned above, steps S105, S106 and S107 can all be performed on the phase difference map, or on the temperature difference map. That is, the step of eliminating the phase drift caused by the magnetic resonance system is performed on the phase difference map and/or the temperature difference map, and the step of correcting the magnetic susceptibility is carried out on the phase difference map and/or the temperature difference map. The step of correcting the phase error caused by motion is performed on the phase difference map and/or the temperature difference map.

Still referring to Figure 3, Figure 3(c) shows the phase difference map for unwrapping and drift correction, Figure 3(d) shows the temperature map, Figure 3(e) shows the susceptibility corrected image, and Figure 3(f) shows the motion error Corrected image.

On the basis of the above embodiments, in another embodiment of the present application, the obtaining the temperature map according to the temperature difference map includes:

S1031: Calculate temperature by using the reference phase difference map and the calibrated phase difference map, and weight the calculated temperature to obtain a temperature map of the target to be measured;

or

A weighted average is performed on the reference phase difference map and the calibrated phase difference map to obtain an average temperature difference, and a temperature map of the target to be measured is calculated according to the average temperature difference.

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 this temperature map is 1, and the weighting coefficients of other phase temperature maps is 0.

After step S1031, it may also include:

S108: Perform multiple interpolation processing on the temperature map of the target to be measured, and use the temperature map of the target to be measured after interpolation processing to calculate an ablation region boundary.

The purpose of performing multiple interpolation processing on the temperature map of the target to be measured is to obtain a smoother boundary of the ablation region, and the specific number of difference processing may be 2 or 3 times.

In the process of calculating the boundary of the ablation region, the following formula is specifically used:

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 τ, and t is the current time. Pixels whose Ω value exceeds a set threshold (eg, 1) are considered ablated.

The ablation calculation method provided in the embodiment of the present application is verified below in combination with specific experiments.

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

As shown in Figure 5, two MR-compatible fiber optic temperature probes were also inserted into the tissue mimic with the probe tips positioned close to the ablation fibers to obtain the gel temperature at various points. Since the fiber optic probe is affected by the ablation fiber during heating, the thermometer only monitors the cooling phase. Figure 5. Specifically, the ablation fiber optic and two fiber optic temperature probes inserted in the gel tissue mimic are shown in Fig. 5, and a gel-filled reference tube is fixed around it as an insulating reference.

Ex vivo experiments with pork and porcine brains were performed using the same scan parameters as the tissue mimic experiments. Two experiments were performed on each type of tissue (gel, pork, pig brain), one with several cycles of laser heating and one with continuous heating and cooling. Calculate the root mean square error between the MR measurement temperature and the fiber measurement temperature as the measurement value of the temperature accuracy.

In vivo experiments in Doberman Pinscher dogs have been approved by the Ethical Review Board of Tsinghua University. Nine adult Doberman pinschers underwent laser interstitial hyperthermia. The heating process was monitored on a 3T MR scanner (Ingenia, Philips Healthcare, Best, The Netherlands) using a multi-echo time gradient echo sequence through 32 receiver coils.

After the ablation procedure, in order to obtain detailed information about the actual extent of the ablation zone, post-ablation images were obtained, including T1 plus gadolinium (T1+Gd) contrast images, fluid-attenuated inversion recovery (FLAIR), diffusion-weighted MR (distortion, by FSL 5.0 by augmentation or corrected by the EPSI method) and T2-weighted images.

Still referring to FIG. 3 and FIG. 4 , FIG. 3 and FIG. 4 illustrate an example of temperature calculation in the ablation calculation method provided by the embodiment of the present application. Fig. 3(a), a time acquired during laser hyperthermia, the coil combination phase image was first obtained by a multi-TE echo sequence; Fig. 3(b), then a phase difference map was obtained, the white arrows indicate the phase around the heating center The phase wrapping that occurs on the map; Figure 3(c) shows the phase difference map after phase unwrapping and B0 drift correction. Fig. 3(d) The temperature map calculated from Fig. 3(c) according to the PRF migration method. White arrows highlight errors due to magnetic susceptibility. Fig. 3(e) Temperature map after magnetic susceptibility correction. White arrows show residual CSF motion-induced errors. Figure 3(f) Temperature map of motion correction.

In Figure 4, an exemplary method flow for a representative pixel includes: Step 1, obtaining a phase difference map and the phase map obtained by unwrapping the reference phase difference map (TE1), as shown by the black arrow, in the case of rapid temperature changes , the phase difference maps corresponding to some echo times are wrapped. Step 2, obtain the phase unwrapping map, step 3, static magnetic field strength (magnetic resonance system) drift correction, that is, B0 drift correction is to reduce system fluctuations, step 4), phase error correction caused by magnetic susceptibility, using the shortest echo Time (TE) corrects for temperature errors caused by susceptibility changes (black arrows) at longer echo times. Step 5 Motion-induced phase error correction. The graphs in the first row and the second row are the phase difference versus time and the corresponding temperature versus time, respectively. The motion-induced phase errors (black arrows) are almost the same for multiple echo times, thus leading to more pronounced temperature errors at shorter TEs. Correcting for motion errors results in smoother phase and temperature profiles.

Tissue mimics and in vitro experimental results:

Figure 6 shows a representative temperature map of isolated pork experiments during laser interstitial hyperthermia. Select six representative images from 300 frames (3 s/frame) acquired during thermal cycling (#50 denotes frame 50, #146 denotes frame 146, and so on). The first row and the second row respectively use the traditional phase unwrapping method and the multi-echo time-based phase unwrapping method (multi-TE unwrapping) proposed in the embodiment of this application. Using state-of-the-art phase unwrapping methods, pixels on the temperature map would be severely damaged due to magnetic susceptibility changes caused by the heat of the laser, irreversible even when the laser was no longer used. The technical principle is as follows: the phase unwrapping method of the prior art is applied to the time dimension for phase jump detection. If the phase difference map of the current frame is unwrapped incorrectly, all subsequent frames will be affected. On the other hand, the ablation calculation method proposed by the present invention is performed on the basis of multi-echo dimensions, thus avoiding interference from previous frames. The third row is a single echo-time-temperature map after phase unwrapping and susceptibility correction, with the damaged pixels around the heating center correctly recovered. The last row is the combination result of multi-echo time data using the ablation calculation method provided by the embodiment of the present application. The resulting MR thermography showed a more uniform temperature at the hot spot.

Figure 7 shows the time-dependent temperature during experiments with gel tissue mimics (Fig. 7(a)) or isolated pork (Fig. 7(b)), measured by two thermometric optical fibers (red lines) and Calculated by the method provided in the embodiment of the present application (dotted black line). In the case of multiple heating (Fig. 7(a)) or single heating (Fig. 7(b)), the temperature-time behavior calculated by the proton resonance frequency (PRF) is compared with that of the fiber optic temperature probe (also known as thermometry fiber optic) measured temperature-time matches very well. Table 2 lists the root mean square error (RMSE) values between the MR calculations and thermometric fiber measurements, which represent the temperature accuracy of the proposed algorithm. Experiment 1 was heated with several laser cycles, while Experiment 2 was a continuous heating and cooling phase. The results showed that the RMSE error for gel, pork or porcine brain tissue was less than 0.5°C in most cases. In Figure 7, the probe (left) represents the temperature measuring fiber on the left, and the probe (right) represents the temperature measuring fiber on the right.

Table 2. Comparison between the temperature measured by the optical fiber and the temperature calculated by MR in the method provided by the embodiment of the present application.

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

Figure 8 shows a representative temperature map of dog 01 in the in vivo experiments. It should be noted that the ablation area is located close to the third ventricle and lateral ventricle. 100 frames (3s/frame) images acquired during laser ablation were selected to be superimposed on T2w magnetic resonance thermal imaging after ablation. From top to bottom is the temperature map calculated from the single echo time (TE) data (TE=6ms and TE=24ms) respectively by the prior art algorithm, and according to the multiple (joint) TE echoes using the algorithm proposed by the present invention. Temperature map computed from the wave sequence. The first row (TE = 6 ms) shows pseudo-hyperthermia in the third lateral ventricle and intraventricular, indicating that short TE calculated temperatures are heavily influenced by CSF flow artifacts. CSF-induced third intraventricular artifacts (indicated by white arrows) still exist in the second row (TE = 24 ms), but are well suppressed by the proposed algorithm for joint TE echo sequences. The second row shows that a longer TE (TE = 24ms) can provide smoother boundaries and better temperature SNR than a shorter TE (TE = 6ms), but as mentioned above, due to the change in magnetic susceptibility , the pixels around the heating center will be damaged. On the other hand, our proposed method integrates the information of multiple echoes, thus the obtained temperature map eliminates both the error caused by CSF and the error caused by magnetic susceptibility, showing a more uniform and symmetrical heating area.

Figure 9 shows the results of laser interstitial thermotherapy ablation in 3 representative dogs (dogs 01-03). The first column is the T2w image after thermal ablation and after administration of gadolinium (post- _T2 ). The location of the ablation fiber can be clearly shown on T2w MRI. The second column shows the final estimated ablation lesions for a given laser dose and the method proposed by the present invention, which are superimposed on the post-ablation T2w in red (darker color on the grayscale image) on the image. The good agreement between the estimated ablation calculated by the method proposed in the present invention and the post-ablation MRI is thus demonstrated. The area of the ablation zone estimated by the method of the present invention is displayed in the upper left corner of the ablation image. The third column is a representative temperature map obtained when the laser heating was most intense. The last three columns are post-ablation FLAIR, DWI and T1w images after gadolinium administration, respectively. They both show a sharp transition zone between dead and living tissue.

Figure 10 shows post-ablation T2w images and final algorithm-predicted damage estimates for all other six cases (dogs 04-09) with different ablation laser doses. Ablation calculations (second row) match well with post-ablation evaluations (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 area of laser ablation ranges from less than 30 mm2 to nearly 90 mm2.

The above experimental results show that the error caused by the magnetic susceptibility in the proton resonance frequency temperature map caused by heating the laser itself can be corrected by using the ablation calculation method provided by the embodiment of the present application. We first propose the application of multi-echo temporal gradient echo pulse sequences to magnetic resonance thermography via the proton resonance frequency shift method. Compared to single-echo sequences, multiple gradient-echo sequences provide more information without requiring additional scan time, and offer new methods for phase unwrapping and artifact removal.

Shorter echo times can tolerate susceptibility artifacts but are sensitive to noise, while longer echo times have better temperature sensitivity and signal-to-noise ratio but are highly affected by susceptibility artifacts. The ablation calculation method proposed by the invention combines the advantages of different echoes to obtain better temperature map measurement results. Moreover, the ablation calculation method of the present invention can significantly improve the robustness and signal-to-noise ratio of magnetic resonance thermography, thereby avoiding damage to healthy tissue caused by misestimating the low temperature.

The method of the present invention also has excellent immunity to CSF flow errors and can provide accurate temperature measurements in or around the ventricle. Compensation for errors caused by CSF 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 computations, making it ideally suited for magnetic resonance thermography where very close to real-time temperature maps are required.

The magnetic resonance temperature imaging system provided in the embodiment of the present application is described below, and the magnetic resonance temperature imaging system described below may be referred to in correspondence with the ablation calculation method described above.

Correspondingly, the embodiment of the present application also provides an ablation calculation system, which includes an ablation calculation module, and the ablation calculation module can execute the ablation calculation method of the present invention. Further, the ablation calculation module may further include:

An information collection module, configured to receive or collect magnetic resonance information, the magnetic resonance information at least includes a phase map corresponding to the echo time obtained by scanning the target to be measured using a gradient echo sequence containing i different echo times, The i is a positive integer greater than or equal to 2;

A temperature difference calculation module, which is used to select at least two groups of phase diagrams corresponding to different echo times to obtain corresponding temperature difference diagrams;

A temperature map calculation module, which is used to obtain a temperature map according to the temperature difference map, and the temperature map may be a weighted temperature map;

An ablation calculation module, configured to calculate the ablation condition of each pixel according to the temperature map.

a phase drift correction module for performing the step of eliminating phase drift caused by the magnetic resonance system,

a susceptibility error correction module for performing the step of correcting susceptibility-induced errors,

A motion error correction module is configured to perform a step of correcting a phase error caused by motion.

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

The embodiment of the present application also provides another magnetic resonance temperature imaging system, including:

A data transmission module, which is configured to receive a magnetic resonance sequence image and judge image integrity;

A temperature calculation module, which is configured to select the sequence, calculate the phase difference, calibrate the phase difference, and calculate the temperature;

A temperature display module, which is configured to display the temperature in the form of a pseudo-color map or an isotherm;

an ablation calculation module, which is configured to calculate and display an ablation result;

Wherein, the time for the system to perform a complete calculation does not exceed 1 second.

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

Correspondingly, the embodiment of the present application also provides a magnetic resonance temperature imaging system, including: a memory and a processor;

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

Correspondingly, an embodiment of the present application further provides a storage medium, on which a program code is stored, and when the program code is executed, the ablation calculation method described in any one of the above-mentioned embodiments is implemented.

To sum up, the embodiment of the present application provides an ablation calculation method and related devices. The ablation calculation method is not a retrospective algorithm or an iterative algorithm, and has a small amount of calculation. It can provide an almost real-time temperature map and has a high D.

The features recorded in the various embodiments in this specification can be replaced or combined with each other. What each embodiment focuses on is the difference from other embodiments. The same and similar parts of the various embodiments can be referred to each other.

The above description of the disclosed embodiments is provided to enable any person skilled in the art to make or use the present 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 the application. Therefore, the present application will not 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 (11)

1. A calculation method for ablation, comprising: Scanning the target to be measured using a gradient echo sequence containing i different echo times to obtain a phase diagram corresponding to the echo time, where i is a positive integer greater than or equal to 2; Selecting at least two groups of phase diagrams corresponding to different echo times to obtain corresponding temperature difference diagrams; Obtaining a temperature map according to the temperature difference map; calculating the ablation condition according to the temperature map; Wherein, obtaining the temperature difference map is carried out as follows: Use the phase diagram at any time to subtract the phase diagram at the reference time to obtain the phase difference diagram at this time, select at least one phase difference diagram corresponding to the echo time as the reference phase difference diagram, and perform other echoes based on the reference phase difference diagram The phase difference diagram to be calibrated corresponding to the time is calibrated to obtain a calibrated phase difference diagram, and the echo time corresponding to the reference phase difference diagram is less than the echo time corresponding to the calibrated phase difference diagram; using the reference phase difference map and the calibrated phase difference map to calculate the temperature difference map at that moment. 2. The ablation calculation method according to claim 1, wherein the calibration is performed as follows: According to the proportional relationship between the phase difference map and the echo time, the estimated value of the phase difference map to be calibrated is calculated based on the echo time and the reference phase difference map; Using the estimated value, unpack the phase difference map to be calibrated according to the phase periodicity to obtain a calibrated phase difference map. 3. The ablation calculation method according to claim 1, further comprising the step of eliminating the phase drift caused by the magnetic resonance system. 4. The ablation calculation method according to claim 3, characterized in that, in the step of eliminating the phase drift caused by the magnetic resonance system, a number of regions with stable physical temperature and uniform tissue are selected as thermal reference points, and are obtained from Phase drift correction is performed by subtracting the average phase difference of the thermal reference point from each phase difference image or subtracting the average temperature change of the thermal reference point from the temperature change map. 5. The ablation calculation method according to claim 1, further comprising the step of correcting the error caused by the magnetic susceptibility, the step of correcting the magnetic susceptibility is carried out on the phase difference map or the temperature difference map, and the step includes: The steps to perform a magnetic susceptibility correction on a temperature difference map include: Obtaining a first temperature map according to the reference phase difference map, obtaining a corresponding second temperature map according to the calibrated phase difference map, judging 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; The temperature value corresponding to the corresponding pixel in the second temperature map is corrected; The steps for susceptibility correction on a phase difference map include: Judging 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 difference map exceeds a preset phase difference threshold, and if so, then The phase difference corresponding to the corresponding pixel in the calibrated phase difference map is corrected. 6. The ablation calculation method according to claim 1, further comprising the step of correcting the phase error caused by motion, the step is by using at least two sets of phase maps corresponding to different echo times at each pixel. Least squares fitting removes motion-induced phase errors. 7. The ablation calculation method according to any one of claims 1 to 6, wherein the temperature map is obtained from a temperature map obtained from at least one reference phase difference map and at least one calibrated phase difference map A weighted temperature map obtained by weighting the temperature map. 8. The ablation calculation method according to claim 1, wherein the ablation is calculated using the following formula: 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 the pixels whose Ω value exceeds the set threshold are considered to be ablation. 9. A storage medium, wherein a program code is stored on the storage medium, and when the program code is executed, the ablation calculation method according to any one of claims 1 to 8 is implemented. 10. An ablation calculation system, characterized by comprising an ablation calculation module capable of executing the ablation calculation method according to any one of claims 1-8. 11. A laser interstitial thermotherapy instrument, characterized in that it comprises: a host, a laser ablation device, and an optical fiber assembly, the host includes a processor, the processor is loaded with program code, and when the program code is executed, it realizes The ablation calculation method according to any one of claims 1-8.