Abstract

Purpose. We present a novel background tissue phase removing method, called anatomical phase extraction (APE), and to investigate the accuracy of temperature estimation and capability of reducing background artifacts compared with the conventional referenceless methods. Methods. Susceptibility variance was acquired by subtracting pretreatment baseline images taken at different locations (nine pretreatment baselines are acquired and called to ). The susceptibility phase data was obtained using the Wiener deconvolution algorithm. The background phase data was isolated by subtracting from the whole phase data. Finally, was subtracted from the whole phase data before applying the referenceless method. As a proof of concept, the proposed APE method was performed on ex vivo pork tenderloin and compared with other two referenceless temperature estimation approaches, including reweighted 1 referenceless (RW-1) and 2 referenceless methods. The proposed APE method was performed with four different baselines combination, namely, (, , , ), (, , , ), (, , , ), and (, , , ), and called APE experiment 1 to 4, respectively. The multibaseline method was used as a standard reference. The mean absolute error (MAE) and two-sample -test analysis in temperature estimation of three regions of interest (ROI) between the multibaseline method and the other three methods, i.e., APE, RW-1, and 2, were calculated and compared. Results. Our results show that the mean temperature errors of the APE method-experiment 1, APE method-experiment 2, APE method-experiment 3, APE method-experiment 4, and RW-1 and 2 referenceless method are 1.02°C, 1.04°C, 1.00°C, 1.00°C, 4.75°C, and 13.65°C, respectively. The MAEs of the RW- and referenceless methods were higher than that of APE method. The APE method showed no significant difference (), compared with the multibaseline method. Conclusion. The present work demonstrates the use of the APE method on referenceless MR thermometry to improve the accuracy of temperature estimation during MRI guided high-intensity focused ultrasound for ablation treatment.

1. Introduction

Recently, high-intensity focused ultrasound (HIFU) has received a surge of attention due to its ability to perform ablation or hyperthermia therapy without an incision [1]. HIFU treatment is superior to conventional resection surgery for several reasons, including reduced cost, shorter recovery time, and greater patient treatment acceptance [2]. During HIFU treatment, monitoring the temperature response is important to ensure that the required thermal dose is delivered to the target region while sparing the surrounding healthy tissues [3]. Magnetic resonance imaging (MRI) is well suited for this purpose because of its great sensitivity of temperature measurement and attractive soft-tissue contrast to clearly reveal protein denaturation and lesion formation [4]. Thus, MRI guided HIFU (MRgHIFU) treatment is widely used in several kinds of diseases, including uterine fibroids [5, 6], brain tumors [7], and essential tremor [8].

MR thermometry can be carried out by several kinds of tissue contrasts, such as proton density [9], T1 relaxation time [10], diffusion coefficient of water molecules [11], and proton resonance frequency shift (PRFS) [12, 13]. Of these, the PRFS method is the most widely used approach because of its linear behavior, ease of measurement, and near tissue-type independence [3, 14]. The PRFS method is based on the chemical shift of water protons [12, 15] and utilizes gradient-echo pulse sequences to acquire MR phase images. With this method, one pretreatment phase image is acquired as a baseline and subsequently subtracted from sequentially acquired phase images during heating to obtain a phase difference and compute the thermal mapping image.

In spite of its success, the PRFS method operates under the assumption that temperature variation is the sole contributor to phase variation. However, in practice, there are multiple sources of phase variations beside temperature variation, such as field drift [16], patient motion [17], heterogeneous fat/water distribution [18, 19], cavitation [20], oxygen concentration changes [21], and ultrasound transducer movement [22]. Among these potential error sources, ultrasound transducer movement is the dominant source of bias during MRgHIFU treatment. These error sources can cause misinterpretation of phase variation as temperature elevation, creating a temperature bias. Correction of this temperature bias is essential to avoid ineffective treatment or unexpected denaturation of adjacent healthy tissues.

Several strategies have been proposed to overcome this issue. A typical method is the multibaseline method, which acquires multiple baseline images to form a lookup table and obtains the corresponding image for baseline subtraction during heating, at the expense of a prolonged scan time [23, 24]. Instead of acquiring multiple baseline images, Rieke et al. proposed a 2 referenceless method which acquires baseline images by applying a two-dimensional polynomial fitting of the nonheated region [25], and Grissom et al. introduced a reweighted 1 regression to this approach to enhance its convenience (hereafter referred to as the RW- referenceless method) [26]. However, tissue-inhomogeneity may cause background artifacts in the temperature maps of both referenceless methods, leading to misjudgment of the focal zone position. These artifacts limit the use of the referenceless methods in clinical applications.

In this study, we introduced an anatomical phase extraction (APE) method to the referenceless MR thermometry whereby the background tissue phase data was isolated and removed before applying polynomial fitting. The APE method was designed to eliminate the detrimental effects of tissue-inhomogeneity. The performance of APE method was evaluated in ex vivo pork tenderloin by using the multibaseline method as standard reference and comparing with the conventional referenceless methods, including 2 and RW-.

2. Materials and Methods

2.1. MRgHIFU System Setup

The MRgHIFU system design and its experimental setup are shown in Figures 1(a) and 1(b). A spherical HIFU-transducer (focal length: 12 cm, aperture radius: 8 cm, and operating frequency: 1.2 MHz) was mounted on an arc structure attached to the MRI patient table. A water bag filled the space between the HIFU transducer and the object to be ablated. The ultrasound power attenuation through the water bag to the ablated object was less than 1.5%. The MR-compatible arc structure can house mechanical mechanisms and associated devices, providing positional control along two degrees of freedom in an in-plane coordinate system. System software was developed with C and JAVA programming language and supported the necessary functions to perform MRgHIFU ablation procedures such as treatment plan, treatment execution with power and positioning controls, workstation-to-MRI scanner communication, target localization, temperature measurement, and ablated tissue necrosis monitoring. Although the fiber optic thermocouple is a common ground truth used to validate MR thermometry, the HIFU transmission would lead to temperature bias of fiber optic thermocouple measurement surrounding the focal region. Therefore, in this study, we did not measure exact temperature changes using the thermocouple device. Since our purpose was to develop a referenceless method which can provide the similar performance as the conventional method, we chose the multibaseline method as the reference standard.

2.2. Ex Vivo Experiments

Ex vivo pork tenderloin was used in this experiment. Nine HIFU ablations with varying power were conducted in different areas of the pork tenderloin (detailed spatial locations and power are given in Table 1) to simulate susceptibility variances caused by device repositioning in the clinical setting. The porcine specimen was obtained 2 hours before the ex vivo MR experiment and kept under the room temperature in the scan room. The positioning of experimental setup and anatomical localization scan took approximately 20 minutes, and the scan time for temperature mapping was approximately 18.5 minutes (1112 seconds). Thus, the total scan time for the whole experiment was approximately 38.5 minutes, and the porcine specimen was kept fresh during the experiment. HIFU energy was transmitted at positions one through nine at seconds, respectively, for a 30-second ablation and 60-second cooling time. The time series of HIFU ablation is also shown in Figure 1(c). For PRFS MR thermometry, the spoiled gradient echo sequence was implemented in a 1.5T MRI scanner (Symphony, Siemens, Erlangen, Germany) with the following parameters: , , , acquisition matrix ,  mm,  mm. The ablation points were preselected, and nine pretreatment baseline images were obtained for the conventional multibaseline thermometry method to acquire standard reference temperature images for comparison (detailed positions are described in Table 1 and relative positions shown in Figure 2(a)). The APE method has to choose four of these baselines to isolate the background tissue phase data. All images were processed and analyzed off-line using MATLAB (R2016b, The MathWorks, Inc, Natick, MA, USA).

2.3. Principle

Let refer to the baseline phase image data when the moving objects (including the HIFU transducer, MRI RF coil, and water bag in our case) are at position , we adopted three assumption for the APE method.

Assumption 1. During point-by-point ablation in MRgHIFU, only two kinds of phase signals exist in the baseline phase data: the background tissue phase data from the fixed object, hereafter called , and the susceptibility phase data from moving objects (including the HIFU transducer, MRI RF coil, and water bag in our case), hereafter called . Thus, the baseline phase image data,, can be expressed as follows:

Assumption 2. When the moving objects are repositioned to different locations, is merely translated in the and directions, and the magnitude of does not change. As for , both magnitude and location of would not be affected by the moving objects repositioning (hence, is also referred to as the fixed phase data).

Assumption 3. The deleterious effects of tissue-inhomogeneity which would lead to background artifacts are contained within .

According to Assumption 3, the background artifacts could be circumvented if can be removed before utilizing the referenceless methods. Thus, the core of the APE method is to isolate and remove from whole baseline phase data . To achieve this aim, we have to acquire first then can be derived from given and in Equation (1).

2.4. The Mathematical Model for the Image Degradation Process

To derive from baseline images, a mathematical model for the image degradation process was used in the APE procedures. In our cases, susceptibility-shift process results in degradation of original (thermal) image. Image degradation process can be expressed as Equation (2) and detailed derivate in the appendix.

where and denote the susceptibility phase data from moving objects when the moving objects are at position and , respectively. The distance between position and can be expressed as and (unit: pixel) along and directions, respectively. is point spread function (PSF) recording the devices shift effect between and (the matrix can also be considered as degrading operation), and the asterisk symbol () indicates a two-dimensional spatial convolution. As and “” are known (detailed in next section), Equation (2) could be cast as a deconvolution problem, and the Wiener deconvolution method can be used to acquire .

2.5. Acquisition of
2.5.1. PSF () Established

is a matrix recording the device shift effect between and . In the case of the transducer translated and (mm) in direction and direction, respectively, definition of is

where the unit of and is pixel.

2.5.2. Acquisition of “” Term

Recall that and refer to the baseline phase image data when the moving objects are at position and , respectively. and can be expressed as

As mentioned previously, is the fixed phase signal, so theoretically should equal and . Thus, “” term can be obtained by

In other words, “” term can be acquired by just subtracting two baseline images.

2.5.3. Using Wiener Deconvolution Method to Acquire

Finally, we perform Wiener deconvolution [27] on Equation (2) to obtain as follows:

where and mean the two-dimensional Fourier transform and inverse Fourier transform, respectively. is a specified constant acting as a free parameter for optimizing the result.

2.6. Phase Unwrapping

All images underwent phase unwrapping prior to the APE method. We used Bruce Spottiswoode’s code from the MATLAB Central File Exchange [28], based on the algorithm proposed by Ghiglia et al. [29].

2.7. APE Procedures

The procedures of the APE method are displayed in a flowchart of Figure 3 and are described below.(1)Pretreatment phase images, ,, , and , were acquired where the moving objects were positioned in four different locations(2)With and , is derived according to Acquisition of described above(3)According to Assumption 1, the background tissue phase caused by the fixed object was isolated by subtracting from (4)To address the noises, was subtracted from , , , and (x, z), respectively, to acquire , , , and . According to Assumption 2, we know that the magnitude of these four susceptibility phase data are equivalent, and there are merely translated in and directions. Thus, these four susceptibility phase data were all shifted to the location and taking the mean as following equation to alleviate the noises(5)Finally, we subtract from to isolate the fixed phase data, called

2.8. Conversion to Temperature

By introducing Assumption 1 that MR phase data acquired during ablation, , is comprised of two kinds of phase signals, and , into the PRFS equation proposed by Ishihara et al. [7], we calculated temperature change as follows:

where is the gyromagnetic ratio,) is the temperature-dependent coefficient for tissue, is the amplitude of the static magnetic field, is the phase data acquired during ablation (comprised of and ), and is the baseline phase data corresponding to . As mentioned previously, is the fixed phase signal, so theoretically should equal . Thus, the equation can be simplified as follows:

To acquire the baseline , we applied a fourth-order RW- referenceless method to to remove phase variance in due to temperature elevation. Then, temperature map can be obtained from Equation (11) with and . Note that we applied a fourth-order RW- referenceless method to to remove phase variance. This step would substantially reduce noise and computational errors generated by the Wiener deconvolution.

2.9. Validation of APE Method

To compare the performance of the APE method with other referenceless methods, we derived the temperature maps by using the APE method, multibaseline method, RW- [26], and referenceless method [25]. We used the code provided by Grissom et al. [22] for the RW- referenceless method and the APE method with fourth-order fitting in our data. Based on Rieke et al.’s recommendation [21], we also used fourth-order polynomial fitting in the referenceless method. Our APE method modified the RW- referenceless method by removing prior to polynomial fitting. We applied four different baseline combination sets to perform APE method. Let baselines which used in APE method as , , , and . The four baseline combinations sets are , (3, 5, 2, 6), (7, 5, 8, 4), and (9, 5, 8, 6) used in APE experiment 1, APE experiment 2, APE experiment 3, and APE experiment 4, respectively. Since our purpose was to develop a referenceless method which can provide the similar performance as the conventional method, we chose the multibaseline method as the reference standard. The accuracy of each MR thermometry map was estimated using the temperature map obtained via the multibaseline method.

To visualize the errors in temperature estimation, we subtracted the multibaseline method temperature map from the temperature maps obtained from the other three methods (APE and RW- and referenceless methods). Furthermore, to evaluate the temperature measurement errors of the three methods compared to the multibaseline method, we drew three rectangular regions-of-interest (ROI) and calculatedthe mean absolute temperature difference (or mean absolute error (MAE)) of each ROI. Two-sample -test was performed to investigate the difference of the mean temperature values of ROIs between standard reference and each method. To examine the reproducibility of APE method, we also calculated the MAEs of each ROI of four APE experiments, respectively. Furthermore, we showed the representative temperature maps derived from these experiments and subtracted the temperature map obtained via experiment 1, experiment 2, experiment 3, experiment 4 from experiment 2, experiment 4, experiment 1, and experiment 3, respectively, to visualize the difference between the temperature maps acquired from these experiments.

3. Results and Discussion

3.1. Results

Susceptibility phase data , which is caused by the moving objects at different positions and derived from the APE method outlined above, is shown in Figure 2(b). The shift, or susceptibility variance, can be visualized by comparing the contours of the central, or reference, position (in red) with the contours of other positions (in blue). Although the variance seems nonsignificant, this susceptibility variance would nevertheless result in noticeable and clinically relevant errors in temperature measurement. Figure 2(c) shows the susceptibility phase data, , and the background tissue phase data, , which is separated from whole phase data . The tissue-inhomogeneity contained in and is obviously observed, whereas the remaining becomes much smoother.

The ROI surrounding the ablation points (ROI1 and ROI2) and the background (ROI3) are shown in Figure 4(a). Each ROI contains 12 pixels. The absolute errors in temperature measurement of the 12 pixels within each ROI were averaged to obtain the MAE of each ROI. A total of 230 frames were captured for each ROI, and the MAEs of these frames for each ROI of each method are plotted in Figures 4(b)–4(g). The mean temperature error was derived by averaging all the MAEs (Table 2). The root mean square errors (RMSEs) between multibaseline method and the other methods of each ROI are also shown in Table 2 to compare the performance of each method. The mean temperature errors of the APE method-experiment 1, APE method-experiment 2, APE method-experiment 3, APE method-experiment 4, and RW-1 and 2 referenceless method are 1.02°C, 1.04°C, 1.00°C, 1.00°C, 4.75°C, and 13.65°C, respectively. The MAEs and RMSEs of the RW- and referenceless method were higher than APE experiments, especially in where the proposed APE method achieved the lowest MAEs and RMSEs of mostly less than 1°C. As compared with multibaseline method, both RW-1 and 2 referenceless methods showed significant difference (). On the contrary, there was no significant difference between multibaseline method and each APE experiment ().

Figures 5(a)–5(d) show the temperature maps derived from the multibaseline, proposed APE, RW-, and referenceless thermometry methods, respectively. Figure 5(e) shows spatial profiles through the heated area of the respective temperature maps. The multibaseline and proposed APE methods demonstrated considerable similarity in both the temperature map and the spatial profile. On the contrary, the method resulted in substantial background artifacts and noticeable discrepancies in both the temperature map and the spatial profile when compared to the multibaseline method. Although the RW- method performed better than the method, background artifacts were still remarkable.

To visualize the error in temperature measurement of each method, Figure 6 shows the representative difference in temperature maps acquired by subtracting the temperature map obtained via each method (APE, RW-, and referenceless) from the multibaseline temperature map. The difference map of APE method exhibited negligible errors. In contrast, both RW- and referenceless methods resulted in noticeable background artifacts in the difference map, with the method demonstrating the largest background artifacts. The background artifacts of the RW- method were still notable, despite being smaller compared to the method.

To validate the reproducibility of APE method, Figure 7 shows the representative temperature map of APE experiments 1 to 4. It is obvious that these temperature maps demonstrate considerable similarity, and the noises in their difference maps are negligible. Note that the color bar range of difference maps in Figure 7 has been narrowed to 1°C to 4°C because the difference values are too small to identify.

3.2. Discussion

The APE method improves upon the referenceless thermometry methods by eliminating background artifacts and temperature measurement errors in the heated region caused by tissue-inhomogeneity. By removing the background tissue phase data before implementing the referenceless methods, we reduced the adverse effects of tissue-inhomogeneity and obtained a temperature map comparable to that of the multibaseline method. Moreover, by requiring only four pretreatment baseline images instead of an entire baseline library required by the multibaseline method, our proposed APE method could potentially reduce scan time and improve patient comfortability.

Theoretically, baselines of different methods should be compared with that of the multibaseline method to determine the accuracy of each method. However, our APE method acquires a different baseline compared to other thermometry methods. Generally, baselines comprise of both and regardless of whether they were acquired via scanning or polynomial fitting. Our strategy removes to resolve the undesirable consequences of tissue-inhomogeneity; thus, our baseline and floating images do not contain Therefore, comparing baselines to determine the accuracy of each thermometry method is not feasible in our case. Instead, to assess the accuracy of each method, we qualitatively compared the temperature maps obtained by each method and quantified the temperature measurement errors by calculating the MAEs.

Repositioning of the moving objects results in shifts in or susceptibility variation. Susceptibility variation in point-by-point ablation (also referred to as sequential sonications) causes errors in temperature measurement in the PRFS method. To our knowledge, and always coexist in phase data and are difficult to distinguish. In this study, we developed a method to separate these two phase sources and remove from whole phase data, thereby visualizing . Consistent with Assumption 2, our results indicate that when the moving objects are repositioned to different locations, the magnitude of does not change; rather, is merely translated in the and directions.

Although both RW- and referenceless methods suffer from negative effects of tissue-inhomogeneity, tissue-inhomogeneity has a smaller influence on the RW- method with regard to the area of the background artifact. The RW-1 method optimizes the referenceless method by automatically excluding the heated pixels and noise from the polynomial fitting procedure. This may have partially alleviated the effects of tissue-inhomogeneity by unintentionally excluding some, but not all, pixels affected by tissue-inhomogeneity. Despite the improvements of the RW- method, background artifacts and temperature measurement errors are still nonnegligible, limiting its clinical applications. Thus, isolation and removal are still necessary.

Theoretically, our APE algorithm can separate and using any two baselines. However, we found that would slightly deform along the shift direction, resulting in small inconsistencies in , which should be fixed. To address this problem and diminish noise, we take the mean of , , , and to acquire (as mentioned in APE Procedures’ step 4) and then subtract from to acquire the fixed phase data .

Our study has some limitations. First, the proposed APE method was developed for two-dimensional susceptibility movement, whereas it is not suitable for susceptibility shifts in three dimensions or cylindrical directions. Thus, it is not viable in three-dimensional or cylindrical positioning devices. Currently, the transducers are typically repositioned along two dimensions (e.g., - directions) in most HIFU ablation treatments, so the proposed APE method is applicable for most of clinical use. In further study, it would be beneficial to improve the proposed APE method to a three-dimensional algorithm for providing more clinical feasibility. Second, the shift of moving objects should be detected when the APE method is performed, but there is not a well-established method to detect these shifts to the best of our knowledge. Thus, it is necessary to develop an approach to record the moving objects shifts in MRgHIFU treatment for strengthening its further clinical use.

This study described an APE method which reduces the effect of tissue-inhomogeneity on referenceless methods and improve the accuracy of PRFS-based temperature measurement. This strategy expands the potential clinical applications of the referenceless methods to beyond homogenous tissues. With our APE method, the referenceless methods can potentially be used in clinical MRgHIFU treatment. Compared to the conventional multibaseline method, the APE method demonstrates comparable accuracy and not only avoids unnecessary risk of treatment interruptions arising from patient movements but also decreases preparation time spent obtaining baseline images. Further studies applying the APE method to in vivo experiments are necessary to validate its feasibility.

4. Conclusion

In this work, we proposed the APE method, a strategy to eliminate the detrimental consequences of tissue-inhomogeneity in referenceless methods. The theoretical aspects of the strategy were described, and an ex vivo experiment was performed for qualitative and quantitative evaluation. Our results indicate that the accuracy of the APE method is comparable with the conventional multibaseline method, implying that the APE method could potentially reduce the total scan time, thereby improving patient comfortability and further reducing the risk of treatment interruption.

Appendix

A. Derivation of Equation (2)

From definition of two-dimensional convolution, the left side of Equation (2) is

From Equation (4), we can see that would be nonzero only if

and

Thus, Equation (A.1) can be rewritten as follows:

According to Equation (4)

Equation (A.4) can be equivalently expressed as

According to Assumption 2, is merely a translation of , and and represent the shift distance along and directions between and . Thus, making the substitution gives

Hence, Equation (2) is proved.

Data Availability

The data used to support the findings of this study are available from the first author ([email protected]) upon request.

Conflicts of Interest

To the best of our knowledge, the named authors have no conflict of interest, financial, or otherwise.

Acknowledgments

This article was subsidized by the Ministry of Science and Technology (MOST 109-2221-E-002-046), the National Health Research Institutes (NHRI-110-BN-PP-06), and the National Taiwan University (NTU), Taiwan.