In this paper, we propose a new Bayesian estimation-based algorithm for two dimensional phase unwrapping of discontinuous phase fields with noisy principal values. The proposed algorithm uses the Markov random field models to build the prior distribution, so the unwrapping problem became equivalent to a minimization of an energy function. Our main contribution in this work is to propose a modification to the classical quadratic potential function, which enforces a global smoothness condition, so that the phase jumps, which result from the phase discontinuities, contribution to the energy are mitigated. This was possible throw weighting the classical quadratic potential function by the probability of occurrence of these jumps which decrease when they increase. Theoretically, this probability follows an hypergeometric distribution which can be approximated as a Gaussian one in order to make easier mathematical manipulations. An analytical expression for the minimization automate was derived.