<p>In the field of digital watermarking, three fundamental and interdependent requirements must be satisfied: robustness, invisibility, and payload capacity. Achieving an optimal trade-off among the three requirements remains a significant challenge. This paper introduces a novel statistical image watermarking method that operates in the nonsubsampled contourlet transform (NSCT)-pseudo Jacobi Fourier moments (PJFMs) magnitude domain, wherein a probability density function (PDF) derived from the Weibull-Burr impounded bivariate distribution (WBIBD) is employed. The proposed statistical watermarking framework consists of two main components: embedding and detection. During the embedding phase, the original image is first decomposed using NSCT, followed by the segmentation of high-frequency subbands into non-overlapping blocks. PJFMs is then computed for NSCT coefficient blocks and digital watermarks are embedded into robust NSCT-PJFMs magnitudes. In the detection phase, the robust local NSCT-PJFMs magnitudes are first modeled using the WBIBD, which accurately captures both the marginal distributions and the strong dependencies of these magnitudes simultaneously. The parameters of WBIBD are computed efficiently by trimmed L-moments estimation. Finally, a watermark detector is then constructed by integrating the WBIBD model with the locally most powerful (LMP) test. Furthermore, closed-form expressions for the watermark detector are derived using the WBIBD framework. Experimental results demonstrate that the probability of detection and AUROC value (area under receiver operating characteristic curve) of the proposed watermarking algorithm are higher than those of other detectors. The results indicate superior detection performance, and it achieves a more effective balance between imperceptibility, robustness, and payload.</p>