We develop a positronium imaging method for the Jagiellonian PET (J-PET) scanners based on the time-of-flight maximum likelihood expectation maximisation (TOF MLEM). The system matrix elements are calculated on-the-fly for the coincidences comprising two annihilation and one de-excitation photons that originate from the ortho-positronium (o-Ps) decay. Using the Geant4 library, a Monte Carlo simulation was conducted for four cylindrical $^{22}$Na sources of β$^{+}$ decay with diverse o-Ps mean lifetimes, placed symmetrically inside the two JPET prototypes. The estimated time differences between the annihilation and the positron emission were aggregated into histograms (one per voxel), updated by the weights of the activities reconstructed by TOF MLEM. The simulations were restricted to include only the o-Ps decays into back-to-back photons, allowing a linear fitting model to be employed for the estimation of the mean lifetime from each histogram built in the log scale. To suppress the noise, the exclusion of voxels with activity below 2% – 10% of the peak was studied. The estimated o-Ps mean lifetimes were consistent with the simulation and distributed quasi-uniformly at high MLEM iterations. The proposed positronium imaging technique can be further upgraded to include various correction factors, as well as be modified according to realistic o-Ps decay models.