Anatomical Phase Extraction (APE) Method: A Novel Method to Correct Detrimental Effects of Tissue-Inhomogeneity in Referenceless MR Thermometry—Preliminary Ex Vivo Investigation

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 φ1 to φ9). The susceptibility phase data φS was obtained using the Wiener deconvolution algorithm. The background phase data φT was isolated by subtracting φS from the whole phase data. Finally, φT 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, (φ1, φ5, φ2, φ4), (φ3, φ5, φ2, φ6), (φ7, φ5, φ8, φ4), and (φ9, φ5, φ8, φ6), 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 t-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- ℓ1 and ℓ2 referenceless methods were higher than that of APE method. The APE method showed no significant difference (p > 0.05), 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.


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 tissuetype 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-ℓ1 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-ℓ1.

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 inplane 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.

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 t = 79, 169, 259, 349, 439, 529, 619, 709, 799 seconds, respectively, for a 30second 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: repetition time ðTRÞ = 37 ms, echo time ðTEÞ = 17:3 ms, flip angle = 18°, acquisition matrix size = 128 × 77, field of view = 256 × 256 mm, slice thickness = 8 mm. The ablation points were preselected, and nine pretreatment baseline images were obtained for the conventional multibaseline thermometry 2 Computational and Mathematical Methods in Medicine 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 offline using MATLAB (R2016b, The MathWorks, Inc, Natick, MA, USA).

Principle.
Let φ j 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 j, we adopted three assumption for the APE method.  Table 1), and thus, the temperature of successive points becomes hotter and hotter.
3 Computational and Mathematical Methods in Medicine data: the background tissue phase data from the fixed object, hereafter called φ T,j , and the susceptibility phase data from moving objects (including the HIFU transducer, MRI RF coil, and water bag in our case), hereafter called φ S,j . Thus, the baseline phase image data, φ j , can be expressed as follows: Assumption 2. When the moving objects are repositioned to different locations, φ S is merely translated in the x and z directions, and the magnitude of φ S does not change. As for φ T , both magnitude and location of φ T would not be affected by the moving objects repositioning (hence, φ T 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 φ T .
According to Assumption 3, the background artifacts could be circumvented if φ T can be removed before utilizing the referenceless methods. Thus, the core of the APE method is to isolate and remove φ T,j from whole baseline phase data φ j . To achieve this aim, we have to acquire φ S,j first then φ T,j can be derived from given φ S,j and φ j in Equation (1).

The Mathematical Model for the Image Degradation
Process. To derive φ S,j 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 φ S,i ðx, zÞ and φ S,c ðx, zÞ denote the susceptibility phase data from moving objects when the moving objects are at position i and c, respectively. The distance between position i and c can be expressed as Δx and Δz (unit: pixel) along x and z directions, respectively. H i−c ðx, zÞ is point spread function (PSF) recording the devices shift effect between φ i and φ c (the matrix H i−c ðx, zÞ can also be considered as degrading operation), and the asterisk symbol ( * ) indicates a two-dimensional spatial convolution. As H i−c ðx, zÞ and "φ S,i ðx, zÞ − φ S,c ðx, zÞ" 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 φ S,i ðx, zÞ.
where the unit of Δx and Δz is pixel.
2.5.2. Acquisition of "φ S,i ðx, zÞ − φ S,c ðx, zÞ" Term. Recall that φ i and φ c refer to the baseline phase image data when the moving objects are at position i and c, respectively. φ i and φ c can be expressed as As mentioned previously, φ T is the fixed phase signal, so In other words, "φ S,i ðx, zÞ − φ S,c ðx, zÞ" term can be acquired by just subtracting two baseline images.

Using Wiener Deconvolution
Method to Acquire φ S,i . Finally, we perform Wiener deconvolution [27] on Equation (2) to obtain φ S,i ðx, zÞ as follows:   Computational and Mathematical Methods in Medicine where FT and FT −1 mean the two-dimensional Fourier transform and inverse Fourier transform, respectively. K 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]. Representative whole phase data (φ 2 ), background tissue phase data (φ T,2 ), and susceptibility phase data (φ S,2 ). The tissue-inhomogeneity of φ T,2 is apparent. When φ T,2 is removed from φ 2 , the remaining φ S,2 is smooth. 5 Computational and Mathematical Methods in Medicine (4) To address the noises, φ T,i ðx, zÞ was subtracted from φ i ðx, zÞ, φ c ðx, zÞ, φ m ðx, zÞ, and φ n (x, z), respectively, to acquire φ S,i ðx, zÞ, φ S,c ðx, zÞ, φ S,m ðx, zÞ, and φ S,n ðx, zÞ. According to Assumption 2, we know that the magnitude of these four susceptibility phase data are equivalent, and there are merely translated in x and z directions. Thus, these four susceptibility phase data were all shifted to the location i and taking the mean as following equation to alleviate the noises (5) Finally, we subtract φ S,iðfinalÞ from φ i to isolate the fixed phase data, called φ T,iðfinalÞ 2.8. Conversion to Temperature. By introducing Assumption 1 that MR phase data acquired during ablation, φ jðabÞ , is comprised of two kinds of phase signals, φ T,jðabÞ and φ S,jðabÞ , into the PRFS equation proposed by Ishihara et al. [7], we calculated temperature change as follows: Wiener filter  Computational and Mathematical Methods in Medicine where γ is the gyromagnetic ratio, α ð−0:01 ppm/°C) is the temperature-dependent coefficient for tissue, B 0 is the amplitude of the static magnetic field, φ jðabÞ is the phase data acquired during ablation (comprised of φ T,jðabÞ and φ S,jðabÞ ), and φ j is the baseline phase data corresponding to φ jðabÞ . As mentioned previously, φ T is the fixed phase signal, so theoretically φ T,jðabÞ should equal φ T,j . Thus, the equation can be simplified as follows: To acquire the baseline φ S,j , we applied a fourth-order RW-ℓ1 referenceless method to φ S,jðabÞ to remove phase variance in φ S,jðabÞ due to temperature elevation. Then, temperature map can be obtained from Equation (11) with φ S,jðabÞ and φ S,j . Note that we applied a fourth-order RW-ℓ1 referen-celess method to φ S,jðabÞ to remove phase variance. This step would substantially reduce noise and computational errors generated by the Wiener deconvolution.

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-ℓ1 [26], and ℓ2 referenceless method [25]. We used the code provided by Grissom et al. [22] for the RW-ℓ1 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 ℓ2 referenceless method. Our APE method modified the RW-ℓ1 referenceless method by removing φ T prior to polynomial fitting. We applied four different baseline combination sets to perform APE method. Let baselines which used in APE method as φ i , φ c , φ m , and φ n . The four baseline combinations sets are ði, m, n, cÞ = ð1, 5, 2, 4Þ, (3, 5, 2, 6), (7,5,8,4), and (9, 5, 8, 6) used in APE experiment  Figure 4: (a) A representative temperature map with the three selected ROIs. Each ROI contains 12 pixels, and the absolute errors in temperature measurement of the 12 pixels within each ROI were averaged to obtain the MAE. (b-g) The MAEs of a total of 230 frames among these three ROIs for the proposed APE #experiment 1 (b), proposed APE #experiment 2 (c), proposed APE #experiment 3 (e), proposed APE #experiment 4 (f), and RW-ℓ1 (d) and ℓ2 referenceless (g) methods, respectively. 7 Computational and Mathematical Methods in Medicine 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-ℓ1 and ℓ2 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 calculated the mean absolute temperature difference (or mean absolute error (MAE)) of each ROI. Two-sample t -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.

Results and Discussion
3.1. Results. Susceptibility phase data φ S , 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 φ S shift, or susceptibility variance, can be visualized by comparing the contours of the central, or reference, position φ S,5 (in red) with the contours of other positions φ S,i (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, φ S,2 , and the Table 2: Mean values and standard deviations for absolute temperature estimation errors, RMSEs, and two-sample t-test analysis between multibaseline and proposed APE, RW − ℓ1 referenceless, and ℓ2 referenceless method, respectively.  Computational and Mathematical Methods in Medicine background tissue phase data, φ T,2 , which is separated from whole phase data φ 2 . The tissue-inhomogeneity contained in φ 2 and φ T,2 is obviously observed, whereas the remaining φ S,2 becomes much smoother. The ROI surrounding the ablation points (ROI 1 and ROI 2 ) and the background (ROI 3 ) 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 methodexperiment 2, APE method-experiment 3, APE methodexperiment 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-ℓ1 and ℓ2 referenceless method were higher than APE experiments, especially in ROI 3 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 (p < 0:05). On the contrary, there was no significant difference between multibaseline method and each APE experiment (p > 0:05).
Figures 5(a)-5(d) show the temperature maps derived from the multibaseline, proposed APE, RW-ℓ1, and ℓ2 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 ℓ2 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-ℓ1 method performed better than the ℓ2 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-ℓ1, and ℓ2 referenceless) from the multibaseline temperature map. The difference map of APE method exhibited negligible errors. In contrast, both RW-ℓ1 and ℓ2 referenceless methods resulted in noticeable background artifacts in the difference map, with the ℓ2 method demonstrating the largest background artifacts. The background artifacts of the RW-ℓ1 method were still notable, despite being smaller compared to the ℓ2 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 φ T 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 φ S and φ T regardless of whether they were acquired via scanning or polynomial fitting. Our strategy removes φ T to resolve the undesirable consequences of tissue-inhomogeneity; thus, our baseline and floating images do not contain φ T : 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 φ S or susceptibility variation. Susceptibility variation in pointby-point ablation (also referred to as sequential sonications) causes errors in temperature measurement in the PRFS method. To our knowledge, φ T and φ S 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 φ T from whole phase data, thereby visualizing φ S . Consistent with Assumption 2, our results indicate that when the moving objects are repositioned to different locations, the magnitude of φ S does not change; rather, φ S is merely translated in the x and z directions.