This paper concerns a basic electrostatic problem: How to calculate generalized Coulomb and self-polarization potentials in heterogeneous dielectric media. In particular, with simulations of ellipsoidal semi-conductor quantum dots and elongated bio-macromolecules being its target applications, this paper extends the so-called three-layer dielectric models for generalized Coulomb and self-polarization potential calculation from the spherical and the spheroidal geometries to the tri-axial ellipsoidal geometry. Compared to the simple step-like dielectric model, these threelayer dielectric models can overcome the mathematical divergence in the self-polarization energy by employing continuous radial dielectric functions. More specifically, in this paper, the novel quasi-harmonic three-layer dielectric model for the ellipsoidal geometry is first discussed, and the explicit analytical series solutions of the corresponding electrostatic problem are obtained in terms of the ellipsoidal harmonics. Then, a robust numerical procedure working for general three-layer dielectric models is developed. The key component of the numerical method is to subdivide the transition layer of the underlying three-layer model into multiple sublayers, and then in each one of them approximate the select dielectric function of the transition layer by one of the quasi-harmonic functional form rather than simply by a constant value as one would normally do. As the result, the numerical method has no numerical divergence.