For any image watermarking algorithm, how to achieve the trade-off among robustness, imperceptibility, and watermark capacity is a challenging problem because of their mutually constrained relationship. In this paper, we design a statistical-based image watermarking system to solve the trade-off problem. In embedding process, consider the imperceptibility and the robustness, we inventively combine the undecimated dual tree complex Wavelet transform (UDTCWT) and the polar harmonic Fourier moments (PHFMs) to obtain the UDTCWT-PHFMs magnitudes as the watermark carriers, and we embed watermark signals in multiplicative manner. In modeling phase, by analyzing the statistical property and the multiple correlations of the UDTCWT-PHFMs magnitudes, we bound the magnitudes of different directions into the vector groups, and then we model them with the vector Roy type bivariate Weibull (BW-Type R) distribution so that we can capture accurately the marginal characteristics, the inter-scale dependencies, and the inter-orientation dependencies at the same time. Moreover, we obtain the model parameters with the double looped maximum likelihood estimation (double-looped MLE). Benefit from the reliable modeling result, we finally derive a novel specific image watermark decoder to blindly extract the hidden watermarks. Extensive experiment results declare the designed statistical image watermarking system achieves the better balance among imperceptibility, robustness, and payload.