PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Magn Reson Med. Author manuscript; available in PMC 2016 April 1.
Published in final edited form as:
PMCID: PMC4224999
NIHMSID: NIHMS587106

Robust Phase Unwrapping for MR Temperature Imaging using a Magnitude-sorted List, Multi-clustering Algorithm

Abstract

Purpose

Several methods in MRI use the phase information of the complex signal and require phase unwrapping (e.g., B0 field mapping, chemical shift imaging, and velocity measurements). In this work, an algorithm was developed focusing on the needs and requirements of MR temperature imaging applications.

Methods

The proposed method performs fully automatic unwrapping using a list of all pixels sorted by magnitude in descending order and creates and merges clusters of unwrapped pixels until the entire image is unwrapped. The algorithm was evaluated using simulated phantom data and in vivo clinical temperature imaging data.

Results

The evaluation of the phantom data demonstrated no errors in regions with signal-to-noise ratios of at least 4.5. For the in vivo data, the algorithm did not fail at an average of more than one pixel for signal-to-noise ratios greater than 6.3. Processing times less than 30 ms per image were achieved by unwrapping pixels inside a region of interest (53 × 53 pixels) used for referenceless MR temperature imaging.

Conclusion

The algorithm has been demonstrated to operate robustly with clinical in vivo data in this study. The processing time for common regions of interest in referenceless MR temperature imaging allows for online updates of temperature maps without noticeable delay.

Keywords: MRI, phase unwrapping, thermometry, multi-clustering, sorted list

Introduction

In MRI, several applications utilize the quantitative phase information of the complex signal and require phase unwrapping, e.g. B0 field mapping, chemical shift imaging, and velocity measurements. The phase unwrapping method, presented in this work, was developed focusing on the needs and requirements of MR temperature imaging (MRTI) applications, which are increasingly used for real-time temperature estimation in monitoring and control of thermal therapy. In particular, several approaches to proton resonance frequency (PRF) shift thermometry, such as referenceless techniques [1,2], require the availability of fast, robust, unwrapped phase maps for temperature estimation.

Numerous algorithms have been proposed to unwrap the phase of MR images. Some methods minimize global cost functions to apply a map of optimal phase offset [3]. These approaches provide high quality results but may result in processing times from several seconds to minutes for images with low signal-to-noise ratios (SNR). Region growing algorithms define unwrapping paths according to quality maps and process images time-efficiently [4,5,6,7]. However, manual parameter adjustments for low SNR images may be required. Fourier techniques unwrap phase fast using a series of transforms [8,9]. Highly wrapped objects on the image boundaries can cause instabilities. Other techniques employ minimum or maximum spanning trees to optimize the unwrapping results [10,11].

Requirements for phase unwrapping algorithms applied in real-time MRTI [5] are the following: 1) The method needs to be fast with respect to the acquisition time to prevent significant delays in temperature monitoring. 2) The algorithm must be robust at low signal-to-noise ratio because during temperature-induced relaxation time changes the MR signal intensity may decrease significantly. Moreover, errors in phase unwrapping may result in major temperature errors. Hence, safe thermal dosimetry might fail. 3) Fully automated algorithms are required because any manual steps during processing are not practical for real-time monitoring.

In this work, a novel fast, robust, fully automated phase unwrapping algorithm appropriate for referenceless PRF thermometry is proposed. The approach provides an alternative method to confine errors in the unwrapped phase maps that might occur near the heat source due to low SNR. Thus, the method is expected to have minimal potential impact on errors in the estimation of the size of the ablation region.

Methods

Background

The measured phase Φmeans(r) ε [0, 2π) of a pixel at a position r ε Rn in an n-dimensional complex MR image I = M · ei · Φmeas, where M is the magnitude, is the wrapped original phase Φorig(r) ε R of the corresponding spin ensemble:

Φmeas(r)=Φorig(r)mod2π
(1)

Phase unwrapping algorithms estimate the offsets 2π · Δ(r), where Δ(r) ε Z is a map of integer values, for reconstruction of the original phase based on the measured phase:

Φorig=Φmeas+2πΔ
(2)

The uncertainty of the phase σ(Φ) is inversely proportional to the signal magnitude M if the signal is sufficiently above background noise [12]:

σ(Φ)=σ(M)M
(3)

The relationship becomes a poor approximation for SNRs around and below 2. However, the phase uncertainty is still monotonically increasing with decreasing signal magnitude even for very low SNRs.

To utilize the relationship between magnitude and phase uncertainty, a sorted list, multi-clustering algorithm is proposed in this work. It performs phase unwrapping using a list of all pixels sorted by magnitude in descending order and creates and merges clusters of unwrapped pixels until the entire image is unwrapped.

Algorithm

Figure 1 shows a flowchart of the proposed phase unwrapping algorithm [13]. First, a list of all m pixels (r1, r2, …, rm) of the complex MR image I is created and sorted according to the corresponding magnitude values M(r1) ≥ M(r2) ≥ (...)M(rm) in descending order. Each pixel ri is processed and corrected as follows.

Figure 1
Flowchart of the phase unwrapping algorithm. First, all pixels are sorted by magnitude to process them in descending order. For pixels without corrected neighbors, new clusters are created. Pixels are added to existing clusters if they are next to one. ...
  • If the pixel ri has no corrected neighbors, a new cluster for this pixel ri will be created and the next pixel ri + 1 in the list will be processed. The weight wi of the new cluster Ci is set to the magnitude M(ri) of the pixel ri.
  • If the pixel ri has corrected neighbors, it will be added to the neighboring cluster Ci,1 with the maximum weight wi,1. The weight wi,1 is the sum of the magnitudes of all pixels inside the cluster.
    wi,1=rjɛCi,1M(rj)
    (4)
    Therefore, the weight wi,1 will be increased by the magnitude M(ri) of the added pixel ri. If the absolute phase difference between pixel ri and its corrected maximum neighbor ri,1 in cluster Ci,1 is more than π:
    |Φmeas(ri)Φmeas(ri,1)|>π
    (5)
    δi,1 ε Z will be selected so that:
    |Φmeas(ri)+2πΔ(ri)Φmeas(ri,1)+2π(Δ(ri,1)+δi,1)|π
    (6)
    This correction step will be completed by updating the value of Δ(ri) to Δ(ri) − δi,1. If there is more than one neighbor with maximum magnitude, Φmeas(ri,1) + 2π · Δ(ri,1) will be replaced by the mean corrected phase of all these neighbors.
  • If the pixel has more neighbors, which are in different clusters Ci,2, Ci,3, …, Ci,k, phase offsets may exist between these clusters. If so, they will have to be corrected. The clusters with weights smaller than the maximum are merged with the cluster with maximum weight. I.e. pixel ri connects existing clusters Ci,1, …, Ci,k of corrected pixels sorted by their weights wi,1wi,2(...)wi,k, where k ε N is the total number of clusters connected in this step. In this case, for each cluster Ci,p, where p ε {2, 3, …, k}, the neighboring pixel ri,p with maximum magnitude is selected. The phase Φmeas(ri,p) of these pixels is used to corrected for phase wraps. If there is more than one neighboring pixel with maximum magnitude in a cluster, Φmeas(ri,p) + 2π · Δ(ri,p) will be replaced by the mean of their corrected phase values. Corresponding offsets δi,2, …, δi,k ε Z are selected, so that:
    ri,p|Φmeas(ri)+2πΔ(ri)Φmeas(ri,p)+2π(Δ(ri,p)+δi,p)|π
    (7)
    The offset values Δ(rq) of all pixels in each neighboring cluster, Ci,p, will be updated to Δ(rq) + δi,p, ∀rq ε Ci,p. The clusters Ci,2, …, Ci,p are merged with cluster Ci,1. The phase of the pixels in cluster Ci,1 with the largest weight wi,1 will not be changed in this step. The next pixel ri+1 then will be processed.

To make the unwrapping paths independent from the positions of the pixels, pixels will be sorted additionally by their phase values if their magnitude values are equal.

In default mode, the algorithm terminates when processing of the entire list of pixels is completed. Additionally, a magnitude threshold can be used to reduce the computation time by stopping early and eliminating the processing of low magnitude pixels (e.g., air), which subsequently can be masked out of the resulting unwrapped phase image.

Implementation

The prototype of this algorithm was implemented in Matlab (R2011; MathWorks, Natick, MA). The code is available at GitHub (https://github.com/fmaier/MRM-2014-PhaseUnwrapping). All processing time measurements were performed using a compiled version of the code on a laptop computer (Dell Latitude E6420, Intel Core i5-2520M, 4 GB RAM). This algorithm processes two-dimensional (2D) and three-dimensional data. However, the algorithm can be easily extended to process n-dimensional data by modifying the neighborhood definition.

Simulated Phantom Data

A simulated ring-shaped phantom (Figure 2) was designed and used to quantify the robustness of the algorithm. The data were stored in a 2D image (256 × 256 pixels). Each sector of the phantom was 30° wide. The signal magnitude decreased from 120 to 10. The phase values were set using Eq. 1 and equation

Figure 2
Simulated phantom with 12 sectors of different SNRs. (a) Magnitude, (b) wrapped phase, (c) unwrapped phase, (d) unwrapping errors. The SNR decreased from 6.0 in sector 1 to 0.5 in sector 12. Close to the center of the phantom the distance between phase ...
Φorig(r)=2rmaxr(r),
(8)

where r(r) is the radius of the pixel from position r to the center of the image and rmax is the radius of the outer boundary of the phantom. Normally distributed noise An external file that holds a picture, illustration, etc.
Object name is nihms587106ig1.jpg(0,400) with a mean value of 0 and a variance of 400 was added to the real and imaginary parts of the complex image. Therefore, the SNR decreased sector by sector in decrements of 0.5 from 6.0 to 0.5.

After application of the algorithm to the simulated data, erroneous pixels were identified by subtracting the original phase and the added noise from the phase unwrapping result. The phase of a pixel was considered to be unwrapped correctly if the absolute difference was smaller than 0.1 π. This simulation was repeated 1000 times, and the mean values and standard deviations (SDs) of the numbers of erroneous pixels per sector were calculated.

For comparison, the PRELUDE algorithm [3] was applied to the same simulated data. The implementation published in the FMRIB Software Library 5.0 (University of Oxford, United Kingdom) was used with default parameters. The algorithm was applied to the entire images.

In vivo Brain Data

The study was conducted with institutional review board approval at The University of Texas M. D. Anderson Cancer Center. PRF temperature imaging [14] datasets from four laser-induced interstitial thermal therapy interventions were used for evaluation of the robustness of the algorithm in in vivo applications. The monitoring images were acquired with an eight-channel, receive-only, intraoperative brain array coil in a 1.5 T open-bore MRI system (Magnetom Espree; Siemens Healthcare, Erlangen, Germany). Figure 3 shows a representative MR temperature image (GRE; flip angle = 20°, TE = 15 ms, TR = 22 ms, matrix = 256 × 128, field of view = 300 × 300 mm2, slice thickness = 5 mm, bandwidth = 260 Hz/pixel). Additionally, the performance of the algorithm on low-SNR images was studied by artificially adding noise to the acquired images. Normally distributed noise was systematically added to the real and imaginary parts of the complex image until the algorithm failed inside the brain region in areas larger than one pixel. The SD was changed in increments of 5.

Figure 3
PRF shift thermometry images of a human head acquired during laser-induced interstitial thermal therapy with different SNRs: (a-e) 22.0, (f-j) 2.7, and (k-o) 2.1. Magnitude images (a, f, and k), wrapped phase images (b, g, and l), unwrapped phase images ...

The apparent SNR was calculated by taking the mean signal magnitude μroi inside the temperature imaging region of interest (ROI) before heating and the standard deviation σair of a 10 × 10 pixels ROI outside the head in air using the following equation [15]:

SNR=0.655μroiσair,
(9)

The in vivo phase unwrapping results were assessed for errors inside the brain region by comparing the unwrapping result with the high-SNR data because the ground truth was not available (compared with the results of the computer simulation on the numerical phantom). The noise simulation was repeated 1000 times for each noise level. No obvious unwrapping errors in the unwrapped high-SNR data were found inside the brain region upon visual inspection.

In total, the algorithm was applied retrospectively to 474 images acquired during thermal ablation procedures to evaluate the processing time. First, the algorithm was applied to the whole image. Second, a magnitude threshold of 5% of the maximum magnitude was employed to prevent unwrapping of noisy air pixels. Third, the algorithm was applied to the ROI (53 × 53 pixels) used for referenceless PRF thermometry. The mean processing time and its SD were calculated.

The PRELUDE algorithm was applied to the in vivo data with the same additional simulated noise. The second processing time measurement was performed using masks based on a magnitude threshold of 5% of the maximum magnitude. In all other experiments, the algorithm was applied to the entire images and ROIs respectively.

Results

Simulated Phantom Data

The results of the application of the algorithm on the simulated phantom data are presented in Table 1. Additionally, Figure 2 shows an example of one experimental run of the algorithm. The first erroneous pixels appear in sector 5 of the phantom at an SNR of 4.0. Phase unwrapping did not fail in sectors 1 – 4, with SNRs of at least 4.5. Large erroneously unwrapped regions were detected in the fourth quarter of the ring phantom at SNRs of no more than 1.5. The low-SNR sectors did not affect the high-SNR sectors. The mean processing time (± SD) required was 0.62 s ± 0.06 s.

Table 1
Total number of erroneous pixels and unwrapping error ratio for each sector of the simulated phantom (Comparison with Figure 2). Mean values and SDs are given for 1000 runs of a random number generator.

In comparison, PRELUDE performed without errors in sectors 1 – 4 as well. In sectors with lower SNR PRELUDE showed fewer errors (Table 1). However, the processing time was 50 s ± 0.6 s.

In vivo Brain Data

For the dataset shown in Figure 3, the measured SNR was 22.0 in the heated region inside the marked ROI. After the addition of normally distributed noise, regions of errors at the boundaries of the brain were found at an SNR of 2.7 (comparison of Figure 3(h-j)). At a lower SNR level of 2.1 large regions of unwrapping errors were found inside the brain (cf. Figure 3(m-o)). Table 2 shows the number of erroneous pixels inside the brain for the tested noise levels. The algorithm failed at an average number of more than one pixel for SNRs lower than 8.1. In comparison, PRELUDE performed similar at SNRs higher than 6.3 (Table 2). At very low SNR, PRELUDE showed fewer errors to the cost of processing time as in the phantom experiment.

Table 2
Total number of erroneous pixels and unwrapping error ratios in the in vivo experiment (Comparison with Figure 3). Mean values, SDs, and processing times are given for 1000 runs for each noise level.

The mean processing time (± SD) required to unwrap entire 2D images was 0.57 s ± 0.03 s. By applying a magnitude threshold of 5% of the maximum magnitude in the image, the algorithm terminated before unwrapping noisy air pixels. Therefore, the processing time was reduced to 0.22 s ± 0.05 s. Faster mean processing times per image of 23 ms ± 3 ms were achieved by unwrapping only pixels in the ROI for referenceless thermometry. In comparison, PRELUDE required 19.1 s ± 3.6 s to unwrap entire 2D images, 0.32 s ± 0.09 s to unwrap masked images, and 134 ms ± 19 ms to unwrap only pixels in the ROI.

Discussion

The sorted list, multi-clustering algorithm allows for fast, robust, fully automatic phase unwrapping and is likely to meet the demands for phase unwrapping for advanced real-time MRTI applications.

The results of the phantom simulation demonstrate good performance of the algorithm with low-SNR data. No errors occurred for SNRs greater than 4. Additionally, because of the order of phase unwrapping, errors in regions with low SNRs are not propagated to regions with high SNRs using this algorithm. Therefore, the algorithm provides very robust results even if the signal magnitude decreases greatly due to heating during thermal therapy. In the case of errors in the center of the heated region, the errors will not propagate to pixels with high SNRs, and valid isotherms at the boundary of the ablation zone still can be created.

Post-processing of the in vivo brain data indicated that the algorithm is also reliable in clinical applications. Errors in the brain region at more than one pixel were found after adding noise to decrease the SNR to 6.3 in the thermometry ROI. The actual SNR of the measured data was about 3.5 times greater than this decreased SNR. Phase unwrapping errors are expected to occur if the spatial phase gradient in the image and the noise add up to an absolute phase difference between two neighboring pixels greater than π. An SNR of 6.3 results in a phase uncertainty (1/SNR) of 0.16 rad. Hence, the uncertainty for the phase differenc (2/SNR) is 0.22. Theoretically, the SNR threshold is 1.8. However, this threshold is valid only if phase differences are caused by noise alone and no other spatial phase gradients exist. Additionally, a large number of pixels might be erroneous if phase unwrapping fails during the merger of two clusters. In areas with very low SNR, this may cause connected regions of erroneous pixels (e.g., sectors 10 – 12 in Figure 2 have such connected regions).

The processing times measured using the compiled Matlab implementation of the algorithm demonstrated the applicability of the algorithm for real-time monitoring of thermal therapy. This phase unwrapping method will cause only minor delays of less than 30 ms per image if the unwrapping is restricted to the thermometry ROI. The increase in time required for unwrapping the phase of the entire MRTI image is mainly caused by noisy air pixels. Also, the increased processing time owed to the creation and merger of many small clusters in this region. However, by using a threshold, the algorithm can be terminated before phase unwrapping of air pixels begins. In future implementations of the algorithm, merging of clusters using efficient data structures such as linked lists should be considered to improve processing times.

The sorted list, multi-clustering algorithm cannot be optimized for parallel computing without major changes, due to the sorted processing of all pixels in descending order. This might limit the use of the algorithm, as large 3D datasets may not be processed in a time efficient way. However, most real-time MRTI algorithms process a limited region in space. Therefore, this algorithm remains a good approach for monitoring thermal therapies.

PRELUDE [3] is a publicly available algorithm and was used as a reference for comparison. The quality of the results was equal for SNRs higher than 4.0 using the simulated phantom data and similar for SNRs higher than 6.3 using the in vivo data. PRELUDE unwraps phase data with very low SNR with fewer errors. However, the processing time increased with noise considerably and was for entire images about 30 – 135 times longer than for the proposed algorithm. ROIs were processed about five times slower. Therefore, the magnitude-sorted list, multi-clustering algorithm is more suitable for real-time monitoring applications.

Compared to classical region growing algorithms [4,5,6,16] this method simultaneously creates and grows regions (i.e., clusters). It always processes the pixel with the least uncertain phase of all remaining pixels. Classical region growing approaches grow regions until checks of quality measures fail. This might result in delayed processing of pixels with high phase certainty due to spatial distance. Additionally, classical region growing approaches may require manual adjustments of the stop conditions depending on the actual image quality.

Images that contain phase singularities caused by some reconstruction algorithms used for parallel imaging might be critical. An irresolvable phase wrap might move in position in an unwrapped image over time, as magnitude values may change for every image in a time series resulting in different pixel processing orders for each image. However, some reconstruction algorithms do not cause such phase singularities in images [17].

Conclusions

The proposed fully automatic magnitude-sorted list, multi-clustering algorithm has been demonstrated to operate robustly with in vivo clinical data. Because of the phase unwrapping order, low-magnitude pixels with noisy phase do not contribute to unwrapping artifacts in high magnitude regions. The processing time for common ROIs in referenceless PRF thermometry of less than 30 ms allows for online updates of temperature maps without noticeable delays.

Acknowledgments

The authors would like to thank Donald R. Norwood (Department of Scientific Publications, The University of Texas M. D. Anderson Cancer Center, Houston, Texas, USA) for editing the manuscript. Valuable discussions on phase unwrapping and its applications with Dr. Wolfgang Stefan, Dr. Renjie He, and Joshua P. Yung (Department of Imaging Physics, The University of Texas M. D. Anderson Cancer Center, Houston, Texas, USA) are also acknowledged.

Grant Sponsor: Portions of this research were supported by NIH grants CA79282, CA016672, 5T32CA119930-03, and 1R21EB010196-01.

References

1. Rieke V, Vigen KK, Sommer G, Daniel BL, Pauly JM, Butts K. Referenceless PRF shift thermometry. Magnetic Resonance in Medicine. 2004;51(6):1223–1231. [PubMed]
2. Salomir R, Viallon M, Kickhefel A, Roland J, Morel DR, Petrusca L, Auboiroux V, Goget T, Terraz S, Becker CD, et al. Reference-Free PRFS MR-Thermometry Using Near-Harmonic 2-D Reconstruction of the Background Phase. IEEE Transactions on Medical Imaging. 2012;31(2):287–301. [PubMed]
3. Jenkinson M. Fast, Automated, N-dimensional Phase-Unwrapping Algorithm. Magnetic Resonance In Medicine. 2003;49(1):193–197. [PubMed]
4. Witoszynskyj S, Rauscher A, Reichenbach JR, Barth M. Phase unwrapping of MR images using ΦUN – A fast and robust region growing algorithm. Medical Image Analysis. 2009;13(2):257–268. [PubMed]
5. de Senneville BD, Maclair G, Ries M, Desbarats P, Quesson B, Moonen CTW. Robust Spatial Phase Unwrapping for On-Line MR-Temperature Monitoring; IEEE International Conference on Image Processing; San Antonio, TX. 2007. pp. III-137–III-140.
6. Cusack R, Papadakis N. New Robust 3-D Phase Unwrapping Algorithms: Application to Magnetic Field Mapping and Undistorting Echoplanar Images. NeuroImage. 2002;16(3):754–764. [PubMed]
7. Ma J. Breath-Hold Water and Fat Imaging Using a Dual-Echo Two-Point Dixon Technique With an Efficient and Robust Phase-Correction Algorithm. Magnetic Resonance in Medicine. 2004;52(2):415–419. [PubMed]
8. Bagher-Ebadian H, Jiang Q, Ewing JR. A modified fourier-based phase unwrapping algorithm with an application to MRI venography. Journal of Magnetic Resonance Imaging. 2008;27(3):649–652. [PMC free article] [PubMed]
9. Schofield MA, Zhu Y. Fast phase unwrapping algorithm for interferometric applications. Optics Letters. 2003;28(14):1194–1196. [PubMed]
10. Ching NH, Rosenfeld D, Braun M. Two-Dimensional Phase Unwrapping using a Minimum Spanning Tree Algorithm. IEEE Transactions on Image Processing. 1992;1(3):355–365. [PubMed]
11. An L, Xiang QS, Chavez S. A fast implementation of the minimum spanning tree method for phase unwrapping. IEEE Transactions on Medical Imaging. 2000;19(8):805–808. [PubMed]
12. Conturo TE, Smith GD. Signal-to-noise in phase angle reconstruction: Dynamic range extension using phase reference offsets. Magnetic Resonance in Medicine. 1990;15(3):420–437. [PubMed]
13. Maier F, Fuentes DTA, Weinberg JS, Hazle JD, Stafford RJ. Proc Intl Soc Mag Reson Med. Salt Lake City, Utah, USA: 2013. Robust Phase Unwrapping Using a Sorted List, Multi-Clustering Algorithm; p. 4294.
14. Rieke V, Butts Pauly K. MR thermometry. Journal of Magnetic Resonance Imaging. 2008;27(2):376–390. [PMC free article] [PubMed]
15. Henkelman RM. Measurement of signal intensities in the presence of noise in MR images. Medical Physics. 1985;12(2):232–233. [PubMed]
16. Yu H, Reeder SB, Shimakawa A, Brittain JH, Pelc NJ. Field map estimation with a region growing scheme for iterative 3-point water-fat decomposition. Magnetic Resonance in Medicine. 2005;54(4):1032–1039. [PubMed]
17. Ros C, Witoszynskyj S, Herrmann KH, Reichenbach JR. Proc Intl Soc Mag Reson Med. Toronto, Ontario, Canada: 2008. Reconstruction of Phase Images for GRAPPA based Susceptibility Weighted Imaging (SWI) p. 1265.