Integrated thermal and magnetic susceptibility modeling for air-motion artifact correction in proton resonance frequency shift thermometry

Abstract Purpose Hyperthermia treatments are successful adjuvants to conventional cancer therapies in which the tumor is sensitized by heating. To monitor and guide the hyperthermia treatment, measuring the tumor and healthy tissue temperature is important. The typical clinical practice heavily relies on intraluminal probe measurements that are uncomfortable for the patient and only provide spatially sparse temperature information. A solution may be offered through recent advances in magnetic resonance thermometry, which allows for three-dimensional internal temperature measurements. However, these measurements are not widely used in the pelvic region due to a low signal-to-noise ratio and presence of image artifacts. Methods To advance the clinical integration of magnetic resonance-guided cancer treatments, we consider the problem of removing air-motion-induced image artifacts. Thereto, we propose a new combined thermal and magnetic susceptibility model-based temperature estimation scheme that uses temperature estimates to improve the removal of air-motion-induced image artifacts. The method is experimentally validated using a dedicated phantom that enables the controlled injection of air-motion artifacts and with in vivo thermometry from a clinical hyperthermia treatment. Results We showed, using probe measurements in a heated phantom, that our method reduced the mean absolute error (MAE) by 58% compared to the state-of-the-art near a moving air volume. Moreover, with in vivo thermometry our method obtained a MAE reduction between 17% and 95% compared to the state-of-the-art. Conclusion We expect that the combined thermal and magnetic susceptibility modeling used in model-based temperature estimation can significantly improve the monitoring in hyperthermia treatments and enable feedback strategies to further improve MR-guided hyperthermia cancer treatments.


Introduction
Hyperthermia is a successful adjuvant to cancer treatments as it improves the effects of chemo-and radiotherapy [1]. During a hyperthermia treatment, the tumor is heated to temperatures ranging between 39 and 44 degrees Celsius for a duration of 60-90 min, which sensitizes the tumor. Noninvasive heating of a tumor can be accomplished by, e.g., Radio Frequency (RF) antennas or Focused Ultrasound (FUS) [2][3][4]. A schematic illustration of a typical treatment setup for RF-based hyperthermia is shown in Figure 1. As the treatment outcome is correlated to the thermal dose in the tumor [5], temperature monitoring and control strategies can improve the treatment quality. Moreover, real-time monitoring of the patient's internal temperature is important given a strong, in-homogeneous, and patient-specific thermoregulatory response.
Monitoring the temperature in clinical practice is typically achieved through intraluminal probes. While these sensors have excellent accuracy and sampling rate, they are spatially sparse and uncomfortable for the patient. To alleviate these issues, Magnetic Resonance Imaging (MRI) can be used to obtain three-dimensional (3D) temperature measurements [6]. However, the so-called Magnetic Resonance Thermometry (MRT) is plagued by a low signal-to-noise ratio (SNR) and low sampling rates [7,8]. To compensate for the low SNR and sample rates, model-based state observers seem a promising solution [9][10][11]. However, possibly the most disrupting measurement corruption, which is currently not corrected for in a model-based observer, is due to the magnetic susceptibility change resulting from airmotion [12,13]. In fact, a recent literature study showed that obtaining an accuracy of 1 degree Celsius is not possible in most regions of the patient [14]. As a result, the limited reliability of MRT hampers its clinical adoption as a replacement for the intraluminal probes despite the benefits of MRT.
Recently, Wu et al. [13] proposed a method for correcting such air-motion-induced susceptibility artifacts. Loosely speaking, [13] employed methods conventionally used in Quantitative Susceptibility Mapping (QSM), namely, the Projection onto Dipole Fields (PDF) [15] and the Laplacian Boundary Value problem (LBV) [16], to remove the dipoleshaped artifacts from thermometry measurements. While both the PDF and LBV algorithm in the method by Wu et al. showed great potential, applying techniques from QSM to thermometry is not straightforward. For example, removing air-motioninduced artifacts using QSM techniques can also result in the removal of part of the temperature-induced phase shift. This issue was partly addressed in [13] by temporarily removing spatially constant and linear terms from the phase measurement prior to removing the air-motion-induced susceptibility artifacts. However, this leaves higher-order terms in the temperature-induced phase shift at risk of being (partially) removed. In particular, the proposed solution is not suitable for scenarios where the heated region is close to air-motion.
To overcome the downsides, in this paper, we propose an integrated thermal and magnetic susceptibility model-based state observer (TSM) to improve the correction of air-motioninduced susceptibility artifacts in MRT. Note that state observers are commonly used to estimate the true state of a system from noisy measurements [17,Ch. 6]. We validate our method using a heated phantom experiment equipped with a movable air volume, see Figure 2. Moreover, we perform an in vivo validation of our TSM method on a data set from a clinical hyperthermia cancer treatment. For both the phantom and in vivo validation, we compare our results to the current state-of-the-art [13].

Theory
In this section, we introduce the thermal and the air-motioninduced susceptibility artifact models which we will use in Section 3.

Thermal model
In this work, we model the thermal dynamics of a patient or phantom using the Pennes' bioheat equation [18], which is a convection-diffusion equation with a distributed source and internal heat loss given by , respectively. Note that the thermal parameters are tissue-specific, possibly even in-homogeneous within a single tissue, and therefore positiondependent. However, the position dependence is omitted in the notation for brevity. Heating is achieved through the absorption of RF waves emitted by m fixed-frequency antennas surrounding the patient, see, e.g., Figure 1. By controlling the power and phase of each antenna the emitted RF waves constructively and destructively interfere and subsequently deposit a localized heat load [19]. Mathematically, the specific absorption rate is given by Here, r [Sm -1 ], E k 2 C 3Â1 [Vm -1 ], p k 2 C, and k Á k 2 denote the electrical conductivity, the complex-valued timeharmonic electric field vector per antenna channel, the complex-valued phasor per antenna channel, and the two-norm, respectively. The time-harmonic electric field generated per antenna for k ¼ 1, . . . , m can be computed using a commercial solver such as Sim4Life [20].
Solving (1) with (2) forward in time is typically performed by spatially discretizing the domain to obtain many coupled ordinary differential equation, which can be solved with wellknown numerical integration techniques, see, e.g., [10]. Throughout this paper we define DTðr, tÞ ¼ Tðr, tÞÀTðr, t 0 Þ as a relative temperature with respect to a baseline at time t 0 , e.g., room temperature for the phantom and approximately 37 degrees Celsius for a patient. In this paper, we consider a homogeneous baseline temperature, however, the presented method is not limited to this particular choice. Note that using relative temperatures is intuitive as the thermometry also measures temperature relative to a baseline.

Air-motion-induced susceptibility artifacts
Hyperthermia cancer treatments are typically monitored using Proton Resonance Frequency Shift (PRFS) thermometry. The PRFS principle exploits the temperature-dependent phase measured by the MRI scanner [7,8]. Besides its   Perfax, which has analogous properties to muscle tissue. In the left cross-section, the tube containing the movable air volume is clearly visible. The black crosses denote four Bowman thermal probes, which are used for validation. temperature sensitivity, the MR signal phase is also dependent on the macroscopic magnetic field [21][22][23], i.e., Duðr, tÞ ¼ uðr, tÞÀuðr, t 0 Þ, (3a) ¼ cTEðB nuc ðr, tÞÀB nuc ðr, t 0 ÞÞ, (3b) %cTEaB 0 ðTðr, tÞÀTðr, t 0 ÞÞ (3c) Note that (4) is based on the assumption that both the macroscopic magnetic field and the susceptibility distribution in (3) do not change over time. However, due to, e.g., magnetic field drift, we need to account for changes in the macroscopic field, denoted by DB mac ðr, tÞ ¼ B mac ðr, tÞÀB mac ðr, t 0 Þ: Additionally, note that the phase shift also depends on the local change in magnetic susceptibility. In this paper, we will neglect the smaller direct influence of the magnetic susceptibility on the phase difference, as we are interested in the more dominant far reaching disturbances that occur with air-motion. To this end, we model the effect of air-motion on the macroscopic magnetic field similar to [22,23] as , and Ã denote the magnetic field drift, the difference in susceptibility distribution, the dipole kernel, a basis vector pointing in the direction of the magnetic field, and the convolution operator, respectively.

Background field removal
In this paper, we will focus on a particular method of QSM to remove the air-motion-induced susceptibility artifacts, namely the PDF method [15]. In this section, we will briefly summarize the underlying principle and assumptions of the PDF method. For a more detailed description, we refer the interested reader to [15] and the references therein.
A key assumption in the PDF method is that the measured temperature field (or equivalently, phase field) has the following decomposition: Here, DTðr, t k Þ [K] denotes the relative temperature difference between the baseline at time t 0 and current time t k , denotes the apparent temperature resulting from B 0 -drift, and DT v ¼ DvÃD a denotes the dipole-shaped apparent temperature resulting from air-motion, which we will refer to as susceptibility artifacts. Note that methods to estimate DT B 0 are well-known [7,8], therefore, we assume that an estimate of DT B 0 is available. Next, the PDF method uses an optimization-based strategy to place susceptibility sources Dv PDF in a pre-defined volume to approximate DT v : Given an approximation of DT v % Dv PDF ÃD a , we can correct the thermometry using Computing the susceptibility sources Dv PDF is performed by solving the optimization problem Here, Dv PDF [-], X S , X ROI denote the identified susceptibility difference distribution, a mask selecting regions with susceptibility artifact sources, and a mask selecting a region of interest in which the clinician would like to measure the temperature, i.e. the patient. The optimization problem (8) can be efficiently solved with the MEDI toolbox [24].

Mask selection
Approximating the background field containing the airmotion susceptibility artifacts using (8) should be done carefully. As mentioned in [15], the PDF method can partly remove the temperature-induced phase shift, i.e. DTðr, t k Þ in (6). Loosely speaking, when X S is chosen too large and includes more than just the susceptibility artifact sources, more of the temperature-induced phase shift is removed. With the same logic, choosing X S too small limits the airmotion artifact removal. Clearly, choosing X S , i.e. the moving air volume, is a crucial aspect in removing air-motion artifacts. To this end, we construct X S by comparing the anatomical scan at the baseline t 0 , and at the current time t k . Any sufficiently large change in intensity between anatomical scans is assumed to indicate air volume movement, and are thus included in X S : Thereto, we compute X S by applying a threshold to the absolute difference between anatomical scans followed by a morphological closing and opening to close holes and remove small objects. More formally, X S is defined as where X 0 S :¼ frjjsðr, t 0 ÞÀsðr, t k Þj ! thresholdg denotes the unprocessed mask, s denotes the signal magnitude corresponding the anatomic MR scan, SðdÞ :¼ frk r k 2 dg denotes a structuring element used to for morphological opening and closing operations, d 1 and d 2 denote the radius of the spherical structuring element, and , § denote the Minkowski sum and difference, respectively. Note that (9) is easily implemented with binary dilation and erosion algorithms found in image processing toolboxes.
The concept of using the difference between two anatomic scans is illustrated in Figure 3 using a phantom equipped with a movable air volume.

Methods
In this section, we start by defining the thermometry acquisition. Hereafter, our TSM air-motion artifact correction scheme is introduced. Next, a heated phantom experiment is detailed, which is used to validate our method against temperature probes. Last, an in vivo validation using data from a hyperthermia cancer treatment is described.

Thermometry acquisition and B 0 -drift correction
The phantom and in vivo temperature is monitored with PRFS-based thermometry acquired by a GE 1.5 T MRI scanner. The MR sequence is a double-echo gradient recalled echo with the first echo at 4.8 ms and the second echo at 19.1 ms. A double-echo scan is used to correct for the electrical conductivity bias [25]. The remaining parameters are: repetition time 620 ms, slice thickness 1 cm, number of slices 25, field of view 50 cm by 50 cm, acquisition matrix 128 Â 128, flip angle 40 , and scan time 82 s. The thermometry acquisition interval is 93 s for the phantom and approximately 10 min in vivo.
The B 0 -drift is compensated by fitting a 3D polynomial, that is linear in all three dimensions, through four oil-filled reference tubes on the exterior of the phantom and through the reference tubes and internal fat in vivo. Both our TSM method and our implementation of [13] will use the same B 0 -drift correction.

TSM air-motion correction scheme
As mentioned in Section 2.3, computing the mask that includes possible susceptibility sources, is a crucial step in the removal of the air-motion artifacts. Our TSM method aims to reduce the removal of the temperature-induced phase shift from the thermometry, when X S is invariably imperfect.
The TSM data processing pipeline is shown in Figure 4. Here, DT measured and DT corrected denote the measured temperature and the air-motion artifact corrected thermometry, respectively. Additionally, DT pred and DT est denote the modelbased temperature prediction and the estimated temperature field that fuses both the thermometry and the modelbased prediction, respectively. To summarize the method in Figure 4, we propose to remove a large part of the temperature-induced phase shift from the thermometry measurement DT measured using model-based temperature prediction DT pred prior to applying the background field removal techniques. Loosely speaking, by subtracting the predicted temperature field from the measurement, there is less temperature-induced phase shift in the measurement for the background field removal techniques to inadvertently remove. After the background field removal, the predicted temperature is added back to obtain the corrected measurement DT corrected : For the remainder of this section, the details of each block in Figure 4 are presented.

B 0 -drift correction
The B 0 -drift in the measured thermometry is corrected by fitting a linear polynomial through reference objects to obtain DT B0 , see Section 3.1 and [7,8].

Prediction
The first step in the model-based observer is the prediction step. Here, the temperature estimate from the previous sample time t kÀ1 is forward simulated using the thermal model to obtain a predicted temperature at time t k . More specifically, we integrate (1) with heat load Q sar ðr, t kÀ1 Þ in (2) over the interval t k Àt kÀ1 starting from the initial state DT est ðr, t kÀ1 Þ: We denote this forward simulation by where f ðÁÞ denotes a numerical integration scheme. In our work, we implement f ðÁÞ based on a continuous-time statespace model (with a finite-dimensional state) obtained by discretizing the spatial domain [10] and subsequently approximate the matrix exponential, which is central in the solution to a set of coupled linear time-invariant differential Figure 3. Illustrating the susceptibility source mask selection X S by detecting air movement using anatomic MR scans. In subfigures A and B, the anatomic baseline scan sðr, t 0 Þ and current scan sðr, t k Þ of the phantom in Figure 2 are shown. Subfigure C shows the absolute difference between A and B. Observe that the moving air is easily distinguished in C. Computing the mask X S is performed by applying an appropriate threshold to C such that all moving air is contained in X S : Figure 4. Schematic overview of the TSM method. The key contribution is the subtraction of the predicted temperature field prior to the background field removal. To obtain the corrected thermometry the predicted temperate is added back to the result of the background field removal.
equations [26]. Note that the particular method of numerical integration is not an intrinsic part of the presented method and alternative integration schemes can be used.

Background field removal
The background field removal is implemented by solving (8) and substituting DT measured ðr, t k Þ by DT measured ðr, t k ÞÀDT pred ðr, t k Þ: After solving (8), the resulting susceptibility distribution is subtracted from the measurement to obtain DT corrected ðr, t k Þ using (7). Note that (7) implicitly accounts for the addition of DT pred as seen in Figure 4. In Figure 5, the intermediate temperature DT measured ðr, t k ÞÀDT B0 ðr, t k ÞÀDT pred ðr, t k Þ, susceptibility field Dv PDF , and identified air-motion-induced artifact DT v are shown to provide better insight into the background field removal method.

Fusion
The last step in our TSM observer is fusing the corrected thermometry with the predicted temperature from the thermal model. A simple policy to combine the predicted and measured temperature is with 0 wðrÞ 1: Intuitively, wðrÞ reflects the trust in either the model-based prediction or the measurement, e.g., wðrÞ%0 in regions with unreliable thermometry and wðrÞ%1 in regions with reliable thermometry. For more details on (11), see Appendix B. Last, note that the estimated temperature DT est ðr, t k Þ will serve as the initial condition for the predicted temperature at time t kþ1 :

Phantom validation
The TSM air-motion artifact correction scheme is validated using a heated phantom experiment. The phantom is equipped with a movable air volume that can be used to inject air-motion artifacts in a controlled manner, see Figure  2 for a schematic illustration. The movable air volume consists of a 10 mm hollow plastic ball attached to a string and inserted inside a tube filled with distilled water. The air volume is moved in between MR scans to introduce air-motioninduced susceptibility artifacts. The resulting estimated temperature field DT est , is validated by comparison to four Bowman temperature probes, as seen in Figure 2. The comparison is quantified using the Mean Absolute Error (MAE) and Standard Deviation (STD) with respect to the probe data for j ¼ 1, . . . , 4: Here, DT probe, j ðt k Þ denotes the relative temperature with respect to time t 0 measured by the j-th probe and r j denotes the position of the j-th probe. Additionally, we compare our method with the current state-of-the-art [13].
Our implementation of the method in [13] will use the same B 0 -drift correction. The background field removal is similar, but instead of subtracting the predicted temperature field, the spatially constant and linear terms are subtracted from the B 0 -drift corrected thermometry. Additionally, the mask X S is chosen as all voxels in the anatomical scan whose magnitude is below a given threshold, as proposed in [13].

In vivo validation
Besides the phantom validation, the TSM method is also validated in vivo using thermometry from a hyperthermia cancer treatment. The patient included in the study is a 43-year-old female diagnosed with primary cervical cancer and treated with curative intent using thermoradiotherapy at Erasmus MC. The patient had a histological confirmed cervical carcinoma with FIGO stage IIIb. The research protocol for this investigation was approved by the medical ethics committee of Erasmus MC, University Medical Center Rotterdam. Note that the approval was obtained prior to the start of the study and the ethics code is MEC 2015-108. Additionally, we have received the informed consent from the patient included in this study.
Similar to the phantom, we placed temperature probes in three intraluminal cavities, namely, the rectum, bladder, and vagina. These probes are used to quantify the performance as detailed in (12). The probe locations with respect to the patients anatomy are shown in Figure 6. Note, that the probes in the patient measure the average temperature along the depth, as seen in right hand side from Figure 6. This behavior is taken into account in the validation. Similar to the phantom validation, we also compare our method to the B 0 -drift corrected thermometry and the method in [13].
The thermal model of the patient is computed using the segmentation in Figure 7. Here, the patient is segmented with four dominant tissues based on MR-scans at the start of Figure 5. Transverse slice of the internal temperature and susceptibility fields used in the TSM background field removal method. Subfigure A shows the difference between the B 0 -drift corrected thermometry and the predicted temperature, i.e. DT measured ðr, t k ÞÀDT B0 ðr, t k ÞÀDT pred ðr, t k Þ: This difference is used to compute the PDF susceptibility distribution. Note that by limiting the PDF method to the difference in measured and predicted temperature, there is less temperature-induced phase shift to inadvertently remove. In subfigure B the identified susceptibility distribution Dv PDF ðr, t k Þ is shown and in subfigure C the identified susceptibility artifact is shown, i.e. DT v ðr, t k Þ ¼ a À1 Dv PDF ðr, t k Þ Ã DðrÞ: Observe that that Dv PDF is confined to a small region defined by X S and note that due to the convolution between the susceptibility distribution and the dipole the impact is observed in larger region than only X S : the treatment [27]. The segmentation used by the thermal model is constant throughout the treatment. Note that a static thermal model can lead to model mismatch as the patients anatomy changes with air motion.

Phantom validation
The probe measurements together with the results of our method, the method in [13], and the B 0 -drift corrected thermometry are shown in Figure 9. A spatial visualization of the thermometry for a single time-step is shown in Figure 8. The performance criteria from (12) applied to the results in Figure 9 are shown in Table 1.
The air-motion-induced susceptibility artifacts are observed in the B 0 -drift corrected thermometry at probe 1 at approximately 30 and 70 min, see Figure 9. Both our TSM method and the method in [13] mitigate the effect of the susceptibility artifact significantly, with improvements in MAE at probe 1 relative to the B 0 -drift corrected thermometry of ten and four times, respectively (Table 1). For the remaining probes, the B 0 -drift corrected thermometry is not corrupted by the air-motion-induced susceptibility artifacts, as seen in Figure 9. We observe that our TSM method is consistent with probes 2, 3, and 4 resulting in a MAE of 0:1 K. The current state-of-the-art yields significantly higher MAE of approximately 0.5 K, as seen by the offset at probes 2, 3, and 4 in Figure 9. Clearly, the method in [13] is susceptible to inadvertently removing the temperature-induced phase shift from the thermometry by the background field removal techniques. Besides the improvements in MAE by the TSM method, we observe a significant reduction in standard deviation, with respect to the B 0 -drift corrected thermometry and the method in [13]. This improvement in standard deviation is the result of the spatial and temporal smoothing from the integrated thermal model.
In addition to the spatially sparse probe comparison, the transverse slice in Figure 8 provides a more visual result of our method. In Figure 8B, observe that the thermometry is corrupted by two dominant susceptibility artifacts. In Figure  8C, the susceptibility artifacts are corrected, while minimally affecting the temperature-induced phase shift. In Figure 8D, the observer-based temperature estimate is shown. Note that the observer can estimate the temperature in regions where no thermometry is available and simultaneously improve the effective SNR by spatially and temporally smoothing the measurement.
For a spatial comparison between our TSM method and the method [13], the temperature along the slice indicated in Figure 8A is shown in Figure 10. From Figure 10, it is clear that the method from [13], partly removed the temperatureinduced phase shift from the thermometry. In contrast, our TSM method is consistent with the thermometry when sufficiently far from the susceptibility artifact ( Figure 10).

In vivo validation
The probe measurements overlaid with the thermometry are shown in Figure 12. Similar to the phantom experiment, a transverse slice for a single time-step is shown in Figure 11. Last, the MAE and STD are provided in Table 2. Observe that the probes are not located in the transverse slice shown in Figure 11, the particular slice location is chosen as it best demonstrated the susceptibility artifact and the corresponding correction.
Notice that, in contrast to the phantom, the in vivo B 0drift corrected thermometry is continuously corrupted with air-motion-induced susceptibility artifacts as, for example, the intestinal gas is continuously in motion, see, e.g., Figures 11  and 12. As a result, a large region of the thermometry is corrupted by susceptibility artifacts, especially in the vicinity of the intestines. Additionally, due to the longer inter-scan time (approximate 10 min) there is less correlation between subsequent scans, compared to the phantom thermometry. As a  result, the temporal smoothing effect of the model-based temperature estimation scheme is reduced.
Based on the quantitative results from Table 2, it is clear that the TSM method outperforms the current state-of-the-art and the B 0 -drift corrected thermometry, both in terms of MAE and STD. Most notably, probes 1 and 2 show a good correspondence with a MAE 0:5 K. The larger discrepancy at probe 3 (MAE %2 K) can be the result of model imperfections in combination with unreliable thermometry. Alternatively, bulk motion or imperfect air-motion correction are also possible explanations. Nevertheless, our method outperforms the B 0 -drift corrected thermometry at probe 3.

Discussion
The phantom with movable air volume provides a controlled experimental setup to verify air-motion artifact correction schemes. This controlled setup clearly showed that the current state-of-the-art potentially removes temperature-induced Figure 8. Transverse slice located at the center of the phantom at 70 min. The subfigures A, B, C, and D denote anatomic scan of the phantom, the B 0 -drift corrected thermometry, the air-motion-corrected thermometry using the TSM method, and the estimated temperature field from the observer. Gray areas indicate that no temperature information is available. Note that in B, several air-motion-induced susceptibility artifacts are present (green arrows) that are removed in C and D. The dashed green arrow in A defines the slice along which temperature is sampled for Figure 10. . Temperature probe correlation for the phantom experiment with B 0drift corrected PRFS thermometry, TSM model-based temperature observer, current state-of-the-art [13]. The light gray area denotes the time-window in which heating was applied. Note that probe 4 is used to calibrate DT measured ÀDT B0 , see Appendix A.  Figure 10. The temperature along the dashed arrow from Figure 8A. The y-axis scale is chosen to demonstrate the severity of the air-motion-induced susceptibility artifact (centered at pixel 95). The gray areas denote regions in which no thermometry is available due to, e.g. air or fat voxels. Observe that the TSM method is successful at removing the susceptibility artifact as the estimated temperature field is consistent with the probes and the B 0 -drift corrected thermometry in regions not corrupted by the susceptibility artifact.
phase shift from the thermometry. This behavior is also suggested by the underlying theory. From Figure 9, it is clear that augmenting the background field removal techniques with temperature field estimates yields a more robust method that is less susceptible to inadvertent removal of part of the temperature-induced phase shift. The performance of our TSM method correlates well with the probe data. Note, however, that this comparison is limited to spatially sparse locations. Nevertheless, we believe that the estimated temperature DT est , using our TSM method is accurate over the whole phantom. The main argument for this, is that the estimated temperature field is consistent with the low spatial frequencies of the thermometry when sufficiently far from susceptibility artifacts, see, e.g., Figure 10 between pixels 20 and 70. This suggest that the background field removal techniques in combination with our model-based state observer prevents removing the temperature-induced phase shift from the thermometry measurements.
The TSM method was also applied to in vivo thermometry from a hyperthermia cancer treatment. Here, we obtained promising results that suggest more of the temperature-induced phase shift is preserved, while simultaneously removing the susceptibility artifacts when compared to the method in [13]. Nevertheless, there are several aspects that could further improve the quality of the estimated temperature field. First, the selection of X S , i.e. the sources of the susceptibility artifacts is performed by comparing the baseline and current anatomical scans. While this works well for the phantom, its application in vivo is less robust. For example, bulk motion could be misidentified as a source for air-motion-induced susceptibility artifacts. Future work on directly detecting the characteristic shape of the susceptibility artifact can improve the robustness of the method. Second, reducing the time in between thermometry scans reduces the impact of limited model quality by increasing the temporal correlation between scans. Third, more sophisticated state observer designs, as seen in [9,10], could result in better temperature estimates from the available measurements. Fourth, more accurate thermal models with, for instance, dynamic segmentation based on the location of the air volume, would further improve the quality of the temperature estimation scheme. Last, in addition to air-motion, in patients with significant bulk movement, we observed significant temperature errors (especially on tissue boundaries), thereby limiting the thermometry quality. Correcting for this type of motion is not included in the TSM model-based temperature estimation scheme. Hence, methods that can compensate for bulk motion can further remove image artifacts from the thermometry. Nevertheless, the in vivo results show that utilizing model-based observers for monitoring hyperthermia treatments is a promising research direction for simultaneous artifact correction and temperature estimation. Figure 11. Transverse slice centrally located in the patient at 60 min. The subfigures B, C, and D denote the B 0 -drift correct thermometry, the air-motion corrected thermometry, and the estimated temperature field from the observer. Gray areas indicate that no temperature information is available. Note that in B, a significant air-motion-induced susceptibility artifact is corrupting a large part of the thermometry (green arrow).

Conclusion
In this paper, we presented an integrated thermal and susceptibility model to effectively remove the air-motion-induced susceptible artifacts in PRFS thermometry using an observer. This TSM method is experimentally validated with a dedicated phantom that allowed for a controlled injection of susceptibility artifacts and in vivo with thermometry data from a hyperthermia cancer treatment. The results obtained in the phantom agreed well with four temperature probes that serve as a ground truth. Additionally, the TSM method substantially improved the standard deviation of the temperature estimate through spatial and temporal smoothing using the integrated thermal model. Similarly, the in vivo results show that the TSM method improves the MAE with respect to the B 0 -drift corrected thermometry and the current state-of-the-art. These results show that, the integrated thermal and magnetic susceptibility modeling using state observers and background field removal methods is a promising research direction that strives to estimate temperature fields from noisy measurements corrupted with motion artifacts.

Disclosure statement
No potential conflict of interest was reported by the author(s).