X-ray emission from the surface of isolated neutron stars (NSs) has been now observed in a variety of sources. The ubiquitous presence of pulsations clearly indicates that thermal photons either come from a limited area, possibly heated by some external mechanism, or from the entire (cooling) surface but with an inhomogeneous temperature distribution. In a NS the thermal map is shaped by the magnetic field topology, since heat flows in the crust mostly along the magnetic field lines. Self-consistent surface thermal maps can hence be produced by simulating the coupled magnetic and thermal evolution of the star. We compute the evolution of the neutron star crust in three dimensions for different initial configurations of the magnetic field and use the ensuing thermal surface maps to derive the spectrum and the pulse profile as seen by an observer at infinity, accounting for general-relativistic effects. In particular, we compare cases with a high degree of symmetry with inherently 3D ones, obtained by adding a quadrupole to the initial dipolar field. Axially symmetric fields result in rather small pulsed fractions ($lesssim 5%$), while more complex configurations produce higher pulsed fractions, up to $sim25%$. We find that the spectral properties of our axisymmetric model are close to those of the bright isolated NS RX~J1856.5-3754 at an evolutionary time comparable with the inferred dynamical age of the source.