Magnetic resonance (MR) phase images reflect various magnetic field variations caused by susceptibility difference, chemical shifts, geometry effects, main field inhomogeneity, etc. [1] Therefore, the phase data has been regarded as very important information to provide fundamental tissue properties complementary to those obtained from the MR magnitude data. However, it has been difficult to quantify and process the phase data because phase wrapping occurs due to the limited dynamic range of phase values by -π to π. There have been conventional phase unwrapping methods based on region-growing methods or high-pass filtering such as branch cut method [2], however, these 2D-based methods often failed to unwrap phase values accurately for regions of complex structures with severe field inhomogeneity or low SNR. Although a voxel-by-voxel unwrapping approach for multi-echo MR images can result in more accurate unwrapped results, it still fails to unwrap some of phase values due to very low-SNR of the multi-echo images. In this study, we propose a robust phase unwrapping method to provide accurate unwrapped phase images for multi-echo MR images based on complex signal modeling.For in vivo experiments, a brain of normal volunteer was scanned by a multi-echo gradient-recalled-echo (MGRE) sequence using a 3T MRI system (Siemens Medical Solutions, Erlangen, Germany). The sequence parameters were the first echo time (TE₁) of 5.67 ms, echo spacing (ES) of 5.51 ms, repetition time (TR) of 95 ms, flip angle (FA) of 27 , slice thickness of 1.6 mm, bandwidth of 444 Hz/Px, field of view (FOV) of 215×215×51.2 mm³, the number of echoes of 16, the number of slices of 32, matrix size of 1024×1024×32×16 interpolated from the data acquired with 512×512×32×16 and in-plane resolution of 0.21×0.21 mm². All types of image processing were performed using MATLAB (The MathWorks, Inc., Natick, MA). The complex signal model used in this study assumes that the phase varies linearly along time. The modeled complex signal $$$\hat{S}\left(\overrightarrow{r},{TE}_{i}\right)$$$ is as follows: $$\hat{S}\left(\overrightarrow{r},{TE}_{i}\right)=\left|{\hat{S}}\left(\overrightarrow{r},{TE}_{i}\right)\right|{e}^{i\left(\omega\left(\overrightarrow{r}\right){TE}_{i}+{\varphi}_{0}\left(\overrightarrow{r}\right)\right)}$$ where $$$\overrightarrow{r}$$$ is a position of a pixel, $$${TE}_{i}$$$ is the ith TE, is the rate of phase change with time of pixel , $$${\varphi}_{0}\left(\overrightarrow{r}\right)$$$ is the phase at TE=0 ms of pixel $$$\overrightarrow{r}$$$, $$$\left|{\hat{S}}\left(\overrightarrow{r},{TE}_{i}\right)\right|$$$ is the magnitude value of $$$\hat{S}\left(\overrightarrow{r},{TE}_{i}\right)$$$. Then by solving the following minimization equation, we can obtain $$$\omega\left(\overrightarrow{r}\right)$$$ and $$${\varphi}_{0}\left(\overrightarrow{r}\right)$$$. $$\min_{\omega\left(\overrightarrow{r} \right),{\varphi }_{0}\left(\overrightarrow{r} \right) }{\sum_{i=1}^{N}{{\left\{{S}_{m}\left(\overrightarrow{r},{TE}_{i}\right)-{\left|{S}_{m}\left(\overrightarrow{r},{TE}_{i}\right) \right| }_{d}{e}^{i\left(\omega\left(\overrightarrow{r} \right){TE}_{i}+{\varphi }_{0}\left(\overrightarrow{r} \right) \right) }\right\} }^{2}} }$$ where $$${\left|{S}_{m}\left(\overrightarrow{r},{TE}_{i}\right)\right|}_{d}$$$ is the denoised magnitude value of pixel $$$\overrightarrow{r}$$$ at $$${TE}_{i}$$$. We used the model-based denoising method to denoise the magnitude data [3]. This denoising method was first applied to the acquired magnitude data . Finally, the unwrapped phase for multi-echo MR images is obtained as follows:$${\varphi}_{u}\left(\overrightarrow{r},{TE}_{i}\right)=\omega\left(\overrightarrow{r}\right){TE}_{i}+{\varphi}_{0}\left(\overrightarrow{r}\right)$$In Fig. 1, we compared the resultant unwrapped phase images for 2D branch cut algorithm, a general thresholding method (thresholding adjacent phase differences by π) and the proposed method. Images in respective three row show phase images from respective different three slices. The phase images with the branch cut algorithm showed many of artifacts due to failure in phase unwrapping. The thresholding method which is a kind of voxel-by-voxel approach resulted in better unwrapped phased images than the bran cut algorithm. In the magnified images (d, j, q), however, many of noise-like scattered artifacts occurred in low-SNR regions, which was due to incorrect unwrapping as seen in black (original) and blue (unwrapped with thresholding) lines in the temporal phase graphs (f, m, t). On the other hand, the phase images with the proposed method do not show any artifacts seen in branch cut or the thresholding method. Also, it can be observed, in the graphs (f, m, t), that the temporal unwrapping were also well performed.We propose a robust phase unwrapping method for low-SNR multi-echo MR images based on complex signal modeling. This method is superior to conventional phase unwrapping methods and provides high-quality unwrapped phase images without any spatial artifacts caused by high noise.