3D Zernike moments based on 3D Zernike polynomials have been successfully applied to the field of voxelized 3D shape retrieval and have attracted more attention in biomedical image processing. As the order of 3D Zernike moments increases, both computational efficiency and numerical accuracy decrease. Due to this phenomenon, a more efficient and stable method for computing high-order 3D Zernike moments was proposed in this study. The proposed recursive formula for computing 3D Zernike radial polynomials combines the recursive calculation of spherical harmonics to develop a voxel-based algorithm for the calculation of 3D Zernike moments. The algorithm was applied to the 3D shape Michelangelo's David with a size of 150×150×150 voxels. As compared to the method without additional acceleration, the proposed method uses a group action of order sixteen orthogonal group and saving unnecessary iterations, the factor of speed-up is 56.783±3.999 when the order of Zernike moments is between 10 and 450. The proposed method also obtained an accurate reconstructed shape with the error rate (normalized mean square error) of 0.00 (4.17×10^-3) when the reconstruction was computed for all moments up to order 450.