<p>Achieving a balance among imperceptibility, robustness, and capacity remains a fundamental challenge in image watermarking. To address this tradeoff, we propose a novel statistical watermarking framework that leverages the nonsubsampled shearlet transform (NSST), the magnitude of the fast polar complex exponential transform (FPCET), and a bivariate Gamma distribution model. We begin by analyzing the robustness and statistical behavior of local NSST-FPCET magnitudes in natural images, revealing notable deviations from Gaussianity and pronounced inter-scale dependencies. These properties are effectively captured by the proposed bivariate Gamma model, which accurately characterizes both the marginal distributions and their joint inter-scale relationships. Model parameters are estimated using restricted maximum likelihood. Building on this statistical foundation, we design a robust watermark decoder based on a maximum likelihood decision rule embedded within the bivariate Gamma framework. Extensive experiments demonstrate that our method outperforms several state-of-the-art statistical and deep learning-based watermarking techniques in terms of robustness and detection accuracy.</p>