This paper addresses interferometric phase (InPhase) image denoising, i.e., the denoising of phase modulo-2p images from sinusoidal 2p-periodic and noisy observations. The wrapping discontinuities present in the InPhase images, which are to be preserved carefully, make InPhase denoising a challenging inverse problem. We propose a novel two-step algorithm to tackle this problem by exploiting the non-local self-similarity of the InPhase images. In the first step, the patches of the phase images are modelled using Mixture of Gaussian (MoG) densities in the complex domain. An Expectation Maximization(EM) algorithm is formulated to learn the parameters of the MoG from the noisy data. The learned MoG is used as a prior for estimating the InPhase images from the noisy images using Minimum Mean Square Error (MMSE) estimation. In the second step, an additional exploitation of non-local self-similarity is done by performing a type of non-local mean filtering. Experiments conducted on simulated and real (MRI and InSAR) datasets show results which are competitive with the state-of-the-art techniques.