A microscopic method for calculating nuclear level density (NLD) is developed, based on the framework of energy density functionals. Intrinsic level densities are computed from single-quasiparticle spectra obtained in a finite-temperature self-consistent mean-field (SCMF) calculation that takes into account nuclear deformation, and is specified by the choice of the energy density functional (EDF) and pairing interaction. The total level density is calculated by convoluting the intrinsic density with the corresponding collective level density, determined by the eigenstates of a five-dimensional quadrupole or quadrupole plus octupole collective Hamiltonian. The parameters of the Hamiltonian (inertia parameters, collective potential) are consistently determined by deformation-constrained SCMF calculations using the same EDF and pairing interaction. The model is applied in the calculation of NLDs of 94,96,98Mo, 106,108Pd, 106,112Cd, 160,162,164Dy, 166Er, and 170,172Yb, in comparison with available data. It is shown that the collective enhancement of the intrinsic level density, consistently computed from the eigenstates of the corresponding collective Hamiltonian, leads to total NLDs that are in very good agreement with data over the entire energy range of measured values.