Quaternion polar complex exponential transform (QPCET) moments and their invariants are widely used as powerful tools in many image processing and pattern recognition applications. However, the accuracy of the conventional approximated method for computing QPCET moments suffers from geometric and numerical errors. This approximated method is very time-consuming. Moreover, computing the high orders of approximated QPCET moments suffers from numerical instability. Computational methods are proposed for fast and accurate computation of the QPCET moments for color images in two coordinate systems. In the first method, the Gaussian quadrature method is applied to compute higher-order moments of QPCET in the Cartesian coordinates. On the other side, an exact kernel-based method is employed to compute the higher-order moments of QPCET in the polar coordinates. A set of numerical experiments is conducted and the obtained results clearly show that the conventional approximated method is unstable, where the numerical instability encountered with moment order , while the first proposed method is unstable for moment order . On the other side, the second proposed method is stable for all orders. The comparison clearly shows the superiority of the second proposed method in terms of image reconstruction capability, numerical stability, fast computation, rotation invariances, and robustness to different kinds of noises.