This paper introduces PnP-DM, a principled Bayesian method for solving imaging inverse problems using diffusion models (DMs) as priors. The method leverages a Markov chain Monte Carlo (MCMC) algorithm based on the Split Gibbs Sampler (SGS) to perform posterior sampling for general inverse problems. By reducing the problem to sampling the posterior of a Gaussian denoising problem, PnP-DM avoids approximations in the generative process and integrates a wide range of state-of-the-art DMs through a unified interface. The method is tested on six inverse problems, including a real-world black hole imaging problem, and demonstrates superior accuracy in reconstruction and posterior estimation compared to existing DM-based methods. Theoretical analysis shows that the average Fisher information distance between the non-stationary and stationary processes decreases with the number of iterations, indicating convergence. Experimental results on both linear and nonlinear inverse problems, including a severely ill-posed black hole imaging problem, confirm the effectiveness of PnP-DM. The method outperforms existing baselines in terms of reconstruction quality and uncertainty quantification. Limitations include computational challenges for large-scale problems and the need for further theoretical analysis of approximation errors. The work has potential applications in computational imaging and related fields, offering a robust solution for image reconstruction with expressive priors.This paper introduces PnP-DM, a principled Bayesian method for solving imaging inverse problems using diffusion models (DMs) as priors. The method leverages a Markov chain Monte Carlo (MCMC) algorithm based on the Split Gibbs Sampler (SGS) to perform posterior sampling for general inverse problems. By reducing the problem to sampling the posterior of a Gaussian denoising problem, PnP-DM avoids approximations in the generative process and integrates a wide range of state-of-the-art DMs through a unified interface. The method is tested on six inverse problems, including a real-world black hole imaging problem, and demonstrates superior accuracy in reconstruction and posterior estimation compared to existing DM-based methods. Theoretical analysis shows that the average Fisher information distance between the non-stationary and stationary processes decreases with the number of iterations, indicating convergence. Experimental results on both linear and nonlinear inverse problems, including a severely ill-posed black hole imaging problem, confirm the effectiveness of PnP-DM. The method outperforms existing baselines in terms of reconstruction quality and uncertainty quantification. Limitations include computational challenges for large-scale problems and the need for further theoretical analysis of approximation errors. The work has potential applications in computational imaging and related fields, offering a robust solution for image reconstruction with expressive priors.