We present a first-principles approach for calculating phonon-polariton dispersion relations. In this approach, phonon-photon interaction is described by quantization of a Hamiltonian that describes harmonic lattice vibrations coupled with the electromagnetic field inside the material. All Hamiltonian parameters are obtained from first-principles calculations, with diagonalization leading to non-interacting polariton quasiparticles. This method naturally includes retardation effects and resolves non-analytical behavior and ambiguities in phonon frequencies at the Brillouin zone center, especially in non-cubic and optically anisotropic materials. It also allows for the natural emergence of hyperbolic polaritons. Furthermore, by incorporating higher-order terms in the Hamiltonian, we also account for quasiparticle interactions and spectral broadening. Specifically, we show how anharmonic effects in phonon polaritons lead to a dielectric response that challenges traditional models. The accuracy and consequences of the approach are demonstrated on GaP and GaN as harmonic test systems and PbTe and $\beta$-Ga${}_2$O${}_3$ as anharmonic test systems.