Phase unwrapping (PU) is essential for various scientific optical applications. This process aims to estimate continuous phase values from acquired wrapped values, which are limited to the interval (-π,π]. However, the PU process can be challenging due to factors such as insufficient sampling, measurement errors, and inadequate equipment calibration, which can introduce excessive noise and unexpected phase discontinuities. This paper presents a robust iterative method based on the plug-and-play (PnP) proximal algorithm to unwrap two-dimensional phase values while simultaneously removing noise at each iteration. Using a least-squares formulation based on local phase differences and reformulating it as a partially differentiable equation, it is possible to employ the fast cosine transform to obtain a closed-form solution for one of the subproblems within the PnP framework. As a result, reliable phase reconstruction can be achieved even in scenarios with extremely high noise levels.