We present a scheme to controllably improve the accuracy of tight-binding Hamiltonian matrices derived by projecting the solutions of plane-wave ab initio calculations on atomic orbital basis sets. By systematically increasing the completeness of the basis set of atomic orbitals, we are able to optimize the quality of the band structure interpolation over wide energy ranges including unoccupied states. This methodology is applied to the case of interlayer and image states, which appear several eV above the Fermi level in materials with large interstitial regions or surfaces such as graphite and graphene. Due to their spatial localization in the empty regions inside or outside of the system, these states have been inaccessible to traditional tight-binding models and even to ab initio calculations with atom-centered basis functions.