On the smallest scales, three-dimensional large-scale structure surveys contain a wealth of cosmological information which cannot be trivially extracted due to the non-linear dynamical evolution of the density field. Lagrangian perturbation theory (LPT) is widely applied to the generation of mock halo catalogs and data analysis. In this work, we compare topological features of the cosmic web such as voids, sheets, filaments and clusters, in the density fields predicted by LPT and full numerical simulation of gravitational large-scale structure formation. We propose a method designed to improve the correspondence between these density fields, in the mildly non-linear regime. We develop a computationally fast and flexible tool for a variety of cosmological applications. Our method is based on a remapping of the approximately-evolved density field, using information extracted from N-body simulations. The remapping procedure consists of replacing the one-point distribution of the density contrast by one which better accounts for the full gravitational dynamics. As a result, we obtain a physically more pertinent density field on a point-by-point basis, while also improving higher-order statistics predicted by LPT. We quantify the approximation error in the power spectrum and in the bispectrum as a function of scale and redshift. Our remapping procedure improves one-, two- and three-point statistics at scales down to 8 Mpc/h.