We present an efficient method for the simulation of time-resolved photoelectron imaging (TRPEI) spectra in polyatomic molecules. Our approach combines trajectory-based molecular dynamics that account for non-adiabatic effects using surface hopping, with an approximate treatment of the photoionization process using Dyson orbitals as initial and Coulomb waves as final electron states. The method has been implemented in the frame of linear response time-dependent density functional theory. As an illustration, we simulate time- and energy-resolved anisotropy maps for the furan molecule and compare them with recent experimental data [T. Fuji, Y.-I. Suzuki, T. Horio, T. Suzuki, R. Mitrić, U. Werner, and V. Bonačić-Koutecký, J. Chem. Phys. 133, 234303 (2010)]. Our method can be generally used for the interpretation of TRPEI experiments allowing to shed light into the fundamental photochemical processes in complex molecules.