A fast referenceless PRFS-based MR thermometry by phase finite difference

Proton resonance frequency shift-based MR thermometry is a promising temperature monitoring approach for thermotherapy but its accuracy is vulnerable to inter-scan motion. Model-based referenceless thermometry has been proposed to address this problem but phase unwrapping is usually needed before the model fitting process. In this paper, a referenceless MR thermometry method using phase finite difference that avoids the time consuming phase unwrapping procedure is proposed. Unlike the previously proposed phase gradient technique, the use of finite difference in the new method reduces the fitting error resulting from the ringing artifacts associated with phase discontinuity in the calculation of the phase gradient image. The new method takes into account the values at the perimeter of the region of interest because of their direct relevance to the extrapolated baseline phase of the region of interest (where temperature increase takes place). In simulation study, in vivo and ex vivo experiments, the new method has a root-mean-square temperature error of 0.35 °C, 1.02 °C and 1.73 °C compared to 0.83 °C, 2.81 °C, and 3.76 °C from the phase gradient method, respectively. The method also demonstrated a slightly higher, albeit small, temperature accuracy than the original referenceless MR thermometry method. The proposed method is computationally efficient (∼0.1 s per image), making it very suitable for the real time temperature monitoring.


Introduction
Magnetic resonance thermometry (MRT) is a promising technique in the guidance and monitoring of many interventional procedures such as tumor ablation (Hynynen and McDannold 2004). Energy delivered to tumors during ablation can be closely monitored by MRT for thermal dose assessment and prevention of tissue charring (McDannold et al 2006). Among the several temperature-sensitive MR parameters previously proposed, the proton resonance frequency shift (PRFS)-based MRT method remains the most popular (Rieke and Butts Pauly 2008) because the PRFS is linear to temperature change and the coefficient is independent of most tissue types (Peters et al 1998).
In the PRFS method, temperature change is calculated from the phase change against a baseline image acquired before heating (Kuroda et al 1997). Motion-induced image misregistration during successive acquisitions is a major source of error for this technique. The multiple baseline method (De Senneville et al 2007) and the model-based referenceless (or self-referenced) method (Rieke et al 2004, McDannold et al 2008 have been proposed to address this problem. In the multiple baseline method, a set of baseline images of the moving object at different positions are acquired before heating. During temperature measurement, the image from the pre-acquired image sets nearest to the imaging slice is selected as the temperature baseline image (Hey et al 2009). This method can also reduce the temperature error in the static region of interest (ROI) caused by the susceptibility distribution change originating from motion nearby (Peters et al 2009). However, it may fail when spontaneous movement occurs and similar baseline image cannot be found in the reference image set. In the original referenceless method (hereby referred as ORG), the baseline phase is found from the image itself by modeling the phase distribution of an unheated region (region of reference, or ROR) as a polynomial using L 2 regression (Rieke et al 2004). The baseline phase of the ROI, where temperature change is to occur, is then extrapolated from the polynomial model (Rieke et al 2004(Rieke et al , 2007. Temperature maps reconstructed from each acquisition independently are therefore insensitive to motion. ORG needs two-dimensional (2D) phase unwrapping before the phase images can be used for polynomial model fitting. Robustness of the phase unwrapping algorithm may be compromised by low SNR or susceptibility artifact (Rieke et al 2004), leading to a large error in the polynomial fitting of phase distribution. Robust 2D phase unwrapping algorithms can be computationally demanding, and would be unsuitable for the time-critical temperature monitoring during thermotherapy. Several referenceless methods have been proposed to address the phase unwrapping issue to improve the computational efficiency. Kuroda et al (2006) modeled the complex MR image with complex polynomials at 0.5 T. Grissom et al (2010a) applied reweighted L 1 regression in the phase finite difference (PFD) image. The zeroth-order term eliminated from PFD of the polynomial is then estimated separately by another L 1 regression-based algorithm. The method does not need manual definition of the heated spot location. Yet, the algorithm requires several user-defined parameters, and its computational load is two orders of magnitude more intensive than the L 2 regression method when fitting the polynomials (Grissom et al 2010a).
Recently an FFT-based phase gradient (PG) method for referenceless MR thermometry that does not need phase unwrapping was proposed (Langley et al 2011). In this PG method, the x-and y-components of the PG images are computed using FFT. The resulting PG images are then fitted to polynomials using L 2 regression. The polynomial model of the baseline phase image is finally derived from the fitted model of the PG image using the relationship described in Liang (1996). However, this FFT-based PG method may encounter two difficulties in a referenceless MRT application. Firstly, the Gibbs ringing caused by the phase discontinuities Figure 1. How the algorithm works. Two vectors along which finite difference and numerical integration of the phase image are to be performed are defined. For a closed ROI boundary, there are two initial conditions for numerical integration. More direction vectors can be used, and they do not need to be orthogonal to each other. at the tissue interface can lead to a large polynomial fitting error. Low pass filtering before the PG calculation by FFT can reduce the ringing effect but it may also introduce a temperature error due to the smoothing effect. Secondly, when recovering phase images from the PG images, the FFT-based method chooses the zeroth-order term (integration constant of the polynomial model) as the average phase in the ROR. When there are regions in the ROR that have total phase wrapping of more than 2π , (not uncommon in high field imaging), estimation of the constant by this method can be far from the truth.
In this paper, a fast referenceless MRT method is proposed to overcome the above two issues in the FFT-based PG method. The proposed method uses PFD to avoid phase unwrapping. The baseline phase image is derived from the PFD image through numerical integration, and the integration constant is derived from the perimeter of the ROI. The performance of the new method will be evaluated by comparing it to PG and ORG through simulation, in vivo room temperature experiments and ex vivo heating experiments.

Referenceless PRFS MR thermometry
In PRFS thermometry, a gradient echo (GRE) sequence is used to acquire the MR temperature change images. The temperature change T is obtained by Kuroda et al (1997) where r is the spatial position of the pixel, ϕ is the phase of the MR image acquired during heating, ϕ b is the baseline phase image acquired before heating, α is the PRFS coefficient, which is −0.01 ppm • C −1 for most biological tissues, γ is the gyromagnetic ratio, B 0 is the magnetic field strength and T E is the echo time of the imaging sequence. In standard PRFS, ϕ b is obtained from a separate scan before any temperature change occurs. In the referenceless method, ϕ b is estimated from the ROR (see figure 1) in the same image by a polynomial model (Rieke et al 2004): (2) Here, the coefficients {C n (m)} are determined from the ROR by L 2 regression. The 2D polynomial is then extrapolated to give the estimate of the baseline phase at the ROI.

Basic principle of the new method
Let I be the complex MR image. The PFD D v (r) of I along a chosen direction v can be defined by where I(r) is the pixel value of the complex image at a spatial coordinate r, ϕ(r) is the phase value at a spatial coordinate r, v is the unit directional vector, conj(.) is the complex conjugate operator and angle(.) is the phase angle operator. The pixel value I(r + v) can be interpolated from the neighboring pixels when it is not on the image grid.
D v (r) in equation (3) gives the PFD image of I along direction v. When the magnitude of the phase difference over adjacent pixels I(r + v) and I(r) is less than π , which is a valid assumption for the temperature change profile, D v (r) computed by equation (3) has no phase wrapping and hence can be used directly for polynomial fitting.
Based on equation (3), the algorithm of PFD can be described as follows. Figure 2 shows the method diagrammatically in a one-dimensional profile for illustration purpose based on a phantom image (figure 2(a)) and a simulated temperature increase (figure 2(b)). Results from PG were also shown step by step for ease of comparison.
(1) Define the polynomial order N, and the several parameters that are required for polynomial model fitting: ROI, ROR and the directional vectors v i (i > 0) (see figure 1 for the definitions of these parameters).
(4) Fit an (N-1)th-order polynomial model (equation (2) above) to the PFD image D i (r) using pixels in the ROR. Note that the polynomial order is reduced by 1, because the finite difference image is a first-order derivative of the phase image to be synthesized. (5) Using the polynomial model in (4), extrapolate the baseline of the PFD image in the ROI, denoted as G i (r). See figure 2(f). (6) Numerically integrate the fitted PFD image G i (r) along v i within the ROI. Note that there are two boundary conditions for a closed boundary in each v i (see figure 2(b)). The resulting baseline phase image θ ij in the ROI related to boundary pixels p ij is given by where p ij is the boundary pixel related to direction v i , j = 1, 2 denotes the two directions of integration for the path v i within the ROI (see figure 1). θ ij (r) is the estimated baseline phase image using the phase of boundary pixel p ij as the initial condition, ϕ(r) is the acquired phase image and C ij is the integration constant that can be found from the boundary condition θ ij (p ij ) = ϕ(p ij ). Figures 2(g) and (h) showed the two phase profiles obtained using boundary conditions A and B, respectively. (7) The temperature induced phase changes in the ROI for a given direction v i and boundary condition p ij can then be computed in the same way as equation (3), as follows: Figures 2(i) and (j) show the results of this operation. (8) The final temperature change image in the ROI is obtained by averaging the phase change images in equation (5): Here, P is the total number of directional vectors and Q is the number of boundary conditions used (typically 2 for computational speed). Note that C ij in equation (4) is the average of the boundary values from all v i used through equation (6) (for instance in figure 1, when two paths are used, four values are averaged). However, averaging is not done immediately after equation (4) in this method. Instead, it is done on ϕ ij (r) obtained from equation (5) instead of θ ij (r) in equation (4). By doing this, any 2nπ phase wrap that may exist in C ij in equation (4) would be canceled through equation (5). The result is that C ij , the constant of integration, is immune to any phase wrap that may come from any part of the image. Note that the baseline phase image is not available from this new method, and that the ROI must lie inside the ROR.
Figure 2(k) shows how the results from the new method compared to the FFT-based PG method in this simulation.

Imaging platform
All experiments were conducted on a 3T MR system (TIM Trio, Siemens, Erlangen, Germany). 2D temperature images were acquired with a GRE sequence using a 1-1 water excitation pulse. Common imaging parameters were: TR = 25 ms, TE =10 ms, slice thickness = 3.0 mm, flip angle = 10 • , bandwidth = 160 Hz/pixel. Image resolution was 1.5 mm × 1.5 mm in all the experiments.

Performance metrics
The temperature error for pixel i in the ROI, denoted as e i , was the difference between the simulated/measured and calculated phase change, converted to degree Celsius ( • C). The error in the whole ROI was characterized using several indices: whereē is the mean error in the whole ROI; (c) absolute mean error AME, which can be derived from (a) and (b) by |ē| = RMSE 2 − SD 2 and (d) 98% percentile of absolute error (PRCT98), which describes the maximal temperature error in the ROI.
These performance indices were used in all the subsequent in vivo and ex vivo heating experiments.

Simulation study
PFD was first compared to PG through a simulation study based on an agar phantom. The phantom image was acquired repeatedly 20 times using the imaging parameters described above. Temperature-induced phase change was created by superimposing a simulated phase change profile (Gaussian distribution, standard deviation = 1.8 pixels, peak value = 0.5 radians and corresponds to a temperature rise of 6.5 • C at 3T with TE = 10 ms) to the acquired phase images. The ROI and ROR were chosen as concentric circles with radii of 18 and 30 pixels, respectively. In order to simulate the inter-scan motion during acquisition, the center of the ROI and ROR was randomly placed on a 10 × 10 grid around the heated area. The simulation was repeated in all the 20 images independently.
The errors in the calculated temperature maps were assessed using the performance metrics in section 3.2.

In vivo experiment
The accuracy of the proposed technique was evaluated in an IRB approved healthy volunteer study performed at room temperature. Six healthy volunteers (all gave informed consent) were recruited. One transverse slice of the brain was imaged (see section 3.1 for imaging parameters). The imaging FOV was 240 mm, the temporal resolution was 3.2 s per slice and the measurement was repeated 40 times.
Since the study took place at room temperature, the calculated temperature maps would exhibit no temperature change. The performance of the new method in the estimation of the baseline temperature image was compared to PG and ORG in terms of RMSE, SD, AME and PRCT98.

Ex vivo heating experiment
The ex vivo heating experiment was performed on bovine liver. Degassed bovine liver was immersed in 1% saline and put in an acrylic container opened at the bottom. The container was fixed in a specially designed platform and put above the high intensity focused ultrasound (HIFU) transducer. Thermoplastic polyurethanes (TPU) membrane was attached to the bottom of the container to separate the bovine liver and saline from the surrounding water (see figure 3). Heating was performed using HIFU as follows. An MR compatible HIFU transducer (focal length/diameter = 10 cm/10 cm at 1.0 MHz, IMASONIC, Besançon, France) was fixed at the bottom of a water tank. The transducer converted the electrical energy delivered from a power amplifier (Amplifier Research, PA, USA) to the ex vivo liver. The applied acoustic power was approximately 25 W in the heating experiments.
Image acquisition was performed during HIFU heating using the parameters in section 3.1 and a FOV of 288 mm. The temporal resolution was 4.8 s per slice. A body matrix and a spine coil were used together for signal reception. The imaging plane was set to the focal plane perpendicular to the ultrasound wave propagation plane (see figure 3).
The experiment consisted of two parts. The first part of the experiment aimed to evaluate the temperature accuracy of the new method compared to PG and ORG in HIFU ablation application. The ex vivo liver was kept stationary. PRFS-based temperature imaging was repeated 60 times. The HIFU transducer was turned on between the 4th and the 40th measurement. This part of the experiment was performed two times (for repeatability).
The second part of the experiment was performed to evaluate the ability of the method in tracking HIFU-induced temperature change during tissue motion. In this part of the experiment, the temperature of ex vivo liver was continuously monitored for 60 measurements. HIFU was turned on at the fourth measurement. At the 14th measurement, the bovine liver was Figure 3. The setup for the ex vivo bovine experiment. An MR compatible HIFU transducer was fixed at the bottom of a water tank. The degassed bovine liver was put in an acrylic container. A piece of TPU membrane separates the bovine liver from water. Imaging plane was set perpendicular to the acoustic wave propagation direction. moved 9 mm along the H-F (head to foot) direction. At the 30th measurement, the bovine liver was moved back to its original position. HIFU heating was stopped during the movement and turned on again before the next measurement. Imaging protocol and coil configuration were the same as the first part of the experiment.
Temperature images from the first part of the experiment were calculated using four methods: the reference subtraction method (or REF), ORG, PG and the proposed method (PFD). REF is liable to the main field drift, especially for long time acquisition. Here, we implemented REF with B 0 drift correction similar to the method used in Salomir et al (2012). The corrected REF method was firstly validated against a fiber optic fluorescent thermometer (Proton Control, Canada) in a water cooling experiment and used as the reference standard. Among the three referenceless methods, the radius of the ROI was chosen to be 14 pixels and the ROR/ROI area ratio was chosen to be 1.0 to ensure that the heating area lies inside the ROI. Note that the center of the ROI was not necessarily the focal point of HIFU. The temperature accuracy of the three referenceless methods compared to REF was assessed using RMSE, SD, AME and PRCT98 for the whole ROI.
Due to motion during HIFU heating in the second part of the heating experiment, comparison was only performed among the three referenceless methods. The ROI-and RORrelated parameters were kept the same as in the first heating experiment.

Data processing
All algorithms were written in MATLAB (Mathworks, NATICK, USA) and run on a computer with 2 GB of RAM and a dual core CPU operating at 2.5 GHz.
Temperature images for both the in vivo and ex vivo experiments were calculated in the following way. For ORG, phase unwrapping before polynomial fitting was done using the quality-guided path finding algorithm (Goldstein et al 1988). During the polynomial fitting procedures in the three referenceless methods, the optimal polynomial order of that specific method, defined as the polynomial order that gives the smallest SD in the first three images of the series over the ROI (Rieke et al 2004), was used. Binary weighting (whereby low signal intensity pixels were masked out) was applied to reduce the effect of noise on the fitting process. In the proposed method, if either I(r) or I(r + v) was masked out, the corresponding pixel in the finite difference image D v (r) was also masked out. Two vectors (v 1 = (1, 0) and v 2 = (0, 1)) were used.

Simulation study
At the simulated hot spot, the temperature calculated from PFD, PG and ORG were 6.36 • C, 6.31 • C and 6.35 • C , respectively, compared to the true temperature increase of 6.5 • C. The accuracies of the three methods assessed in the whole ROI using the four performance indices described earlier were tabulated in table 1. The new method shows improved accuracy compared to ORG and PG. Simulation results showed that the estimated temperature at or around the hot spots from PG and PFD are close to the true values (see figure 2(k)). High temperature errors occurred mainly at the borders of the ROI, and PFD out-performed PG in this region. The average computational time was less than 0.1 s per image for both PG and PFD.

In vivo experiment
The performance of the three referenceless methods based on the six volunteers from the in vivo brain experiment was listed in table 1. Consistent with the simulation results, the proposed method shows obvious accuracy improvement compared to PG. It also shows slightly better accuracy than ORG in every volunteer study.

Ex vivo heating experiment
Over the entire HIFU heating and cooling session (60 measurement), the RMSE/SD/AME/PRCT98 of the proposed method compared to REF from the two ablations were listed in table 1. The new method has the highest temperature accuracy and agrees with the findings from simulation and the in vivo experiment. Figure 4 shows the 2D temperature images from REF and the three referenceless methods at the 40th measurement where the highest temperature rise occurred, overlaid on the anatomical image. At the lower left corner of the ROI, all three referenceless methods showed obvious temperature deviation compared to REF. This may be attributed to the susceptibility change at the tissue-vessel interface around the blood vessel in that region (see arrow). The results suggest that the effect of susceptibility artifact on the three referenceless methods is practically similar. Figure 5 shows the temperature change at the focal point of HIFU heating from the first part of the ex vivo bovine liver heating experiment (bovine liver #1 in table 1). Temperature changes over time at the HIFU focus measured by the four methods are very similar visually ( figure 5(b)). Figure 5(c) shows the temperature profile at the focal point along the white dotted line in figure 4(c) at the 40th measurement when HIFU was just turned off.  Results from the moving table heating experiment are shown in figure 6. The time courses of temperature change at the focal point during the entire HIFU heating process from the three different referenceless methods and REF are shown in figure 6(a). Figure 6(b) shows the temperature difference between the three referenceless methods and REF. The increased temperature differences between the 14th measurement (67 s) and the 30th measurement (144 s) corresponded to the time when the bovine liver was displaced by 9 mm from its original position (and REF is no longer a good reference standard). All the referenceless methods studied were able to track the temperature changes despite the motion. Figure 6(c) shows the temperature maps at three instances of the ablation experiments (marked I, II and III in figure 6(a)).

Discussion
The simulation and the several experiments showed that the new method is consistently more accurate than PG (see table 1), and slightly better than ORG. The improvement of temperature accuracy may be attributed to the use of (a) finite difference instead of FFT for PG calculation and (b) boundary values at the ROI to determine the constant term in the polynomial model eliminated during the calculation of the finite difference images.
While the use of either phase derivative or finite difference images in polynomial fitting can address the phase unwrapping problem in phase images, use of the FFT relation to derive the phase derivative image in PG introduces Gibbs ringing whenever phase discontinuity appears, compromising temperature accuracy in the referenceless MR thermometry application. By using finite difference instead, abrupt phase changes are localized. These points become outliers in the model fitting process and have little impact on the polynomial model used in PFD. This may also explain why the optimal polynomial order was around 2 to 3 in the new method. In fact, low order polynomials for reference MR thermometry are advantageous from the computational viewpoint as they are more robust toward finite word length effects in computers .
The inclusion of boundary conditions in the polynomial fitting process is trivial in PFD (equation (4)). This consideration leads to a more accurate estimation of the integration constant and helps reduce the temperature error compared to PG. It may also contribute to the obvious albeit small reduction in the temperature error of the new method compared to that of the ORG found in this study. In fact, equation (6) may be modified to support weighted averaging, i.e., Our preliminary study found that the above method further improves temperature accuracy in a small but noticeable way.
Magnitude-squared weighting for polynomial fitting is common in L 2 regression algorithms while binary weighting (whereby only magnitude signal higher than a threshold are included for fitting) is more suitable in a low coil sensitivity uniformity situation (Rieke et al 2004). In the in vivo room temperature experiment using a 12-channel head coil which has homogeneous signal sensitivity over the imaging region, results obtained using magnitudesquared weighting were consistent with the binary weighting scheme. However, in the heating experiment, the non-uniform coil sensitivity distribution of the combined body matrix coils and spine coils was suboptimal for magnitude weighting in the polynomial fitting process. As such, binary weighting was used in all experiments here to facilitate comparison of results.
In the new method, better accuracy may be achieved by increasing the number of direction vectors and therefore incorporating more information about the ROI perimeter in equation (6). Two orthogonal vectors (along the readout and phase encoding direction of the MR images) were used in this study for simplicity and computational speed since each vector requires a separate fitting procedure. It was also found that the direction vector's orientation could be arbitrarily chosen. This flexibility is particularly important in practice when part of the ROI boundary may be corrupted by susceptibility artifacts. The PG approach may not be able to provide such flexibility in a trivial way.
Table 1 also shows that the performance metrics decrease in the room temperature experiments and decreases further in the HIFU heating experiment for all the three methods. In the room temperature experiments, the polynomial model was fitted to phase images that contain phase variation coming from tissue susceptibility, flow, etc. In the HIFU experiments, temperature-induced tissue property changes further made the extrapolation of the reference phase map less accurate than the room temperature experiment.
In figure 4, the temperature maps from the referenceless methods were slightly different from that of the reference method, especially in the lower left part of the focal point. Such difference was also found in figure 6(c). Such temperature error may be attributed to the susceptibility change at the tissue-vessel interface around the blood vessel in that region. In referenceless MR thermometry, the phase model built from the ROR does not include the local field disturbance caused by the susceptibility change inside the ROI.
In fact, the temperature-induced susceptibility change is known to be an error source for the PRFS method. For instance, the temperature dependence of fat susceptibility is in the same order as the water proton screening factor (Sprinkhuizen et al 2010). In referenceless MR

Phase Image
Unwrapped Phase Image thermometry, it has been shown that this error may be reduced by excluding the susceptibilityinduced field change from the phase image during model fitting (Kickhefel et al 2012). Hybrid multiple baseline and referenceless methods (Grissom et al 2010b, De Senneville et al 2010 have also been proposed to address this problem and our proposed method does not conflict with these hybrid schemes. Our ex vivo heating experiment showed that the three referenceless methods were immune to inter-scan motion during successive acquisition. For in vivo temperature monitoring, intrascan motion due to respiration and/or heart beat would bring in undesirable phase change or artifact in the temperature maps and reduce the accuracy of the temperature maps. This problem could be addressed by using fast imaging techniques such as echo planar imaging or parallel imaging (Stafford et al 2004). The effect of these fast imaging techniques on temperature accuracy (due to reduced SNR) needs to be carefully considered, and is beyond the scope of this study.
2D phase unwrapping is time consuming. For instance, the Goldstein quality-guided path finding algorithm (Goldstein et al 1988) used in this study took about 5 s for a single 192 × 192 image. The reweighted L 1 regression method (Grissom et al 2010a) does not require phase unwrapping, but the model fitting process is computationally intensive. Computational power may not seem to be an issue nowadays when GPU are available. However, in real time monitoring of thermal therapy where referenceless MR thermometry is only one step among a set of post-processing algorithms (Hey et al 2009), a fast and robust MR thermometry method would allow for more efficient allocation of computational resources to other processes. Results from experiments performed here found that the newly proposed PFD method is an improvement of the FFT-based PG method over temperature accuracy without additional demand on computational power.
On the other hand, the 2D phase unwrapping algorithm may not be very robust. It may fail for a variety of reasons such as low SNR, severe susceptibility artifact or the presence of poles. The results from the in vivo volunteer experiment showed that PFD worked equally well in modeling the complex in vivo phase images in real-life situations. In one volunteer experiment, (not included in the analysis) the Goldstein quality-guided path finding algorithm (Goldstein et al 1988) failed to unwrap the brain image due to the presence of a pole in the phase map (figure 7, black solid arrow). The RMSE/SD/PRCT98 of PFD compared to REF in this case is 1.09/1.09/3.38 • C using the original phase image ( figure 7(a)). However, the RMSE/SD/PRCT98 of ORG compared to REF is 2.28/2.13/5.65 • C using the unwrapped phase image ( figure 7(b)). Both methods used identical ROR/ROI configuration in figure 7(b). In this particular case, the error of the ORG method increases significantly (∼100%) due to the failure of this phase unwrapping algorithm, while the proposed PFD method achieved similar accuracy comparing to the results from the six volunteers (see table 1). In real-life HIFU ablation procedures where the existence of poles may not be known beforehand, the need to avoid phase unwrapping in referenceless MR thermometry becomes more apparent.
While the reweighted L 1 regression algorithm (Grissom et al 2010a) also used the finite difference approach to avoid phase unwrapping, it differs from the proposed method in two important ways. Firstly, they used a method similar to Langley et al (2011) to find the constant of integration when recovering the extrapolated baseline phase images. When the phase difference inside the ROI between the fitted baseline phase without dc and phase image is near −π /+π , the averaged constant value will have large error. In PFD, the constant value is not explicitly computed (see section 2 above), and averaging takes place after the temperature change maps are computed independently according to different boundary pixels. Therefore, the phase wraps are dealt with through equation (5) before they are averaged. Secondly, the reweighted L 1 regression algorithm is used for polynomial fitting to eliminate the need for ROI definition in Grissom et al (2010a) (only the ROR is needed). The procedure requires four user-defined parameters, and the regression procedure is more time-demanding than the L 2 regression method. A comprehensive comparison on referenceless methods considering accuracy, computational requirements and robustness to a low SNR and susceptibility artifact should be explored further in the future.

Conclusion
An improved referenceless method that does not need two-dimensional phase unwrapping has been proposed. Simulation and experiments show that this method is more accurate than the phase gradient method and the original referenceless method with reduced computational requirements. The method will be useful in real time referenceless MR thermometry.