A notorious problem in observational astrophysics is the recovery of the true shape and spectral energy distribution (SED) of a galaxy despite absorption by interstellar dust embedded in the same galaxy. It has been solved only for a few hundred nearby galaxies with exquisite data coverage, but these techniques will not be applicable to the billions of galaxies in upcoming large wide-field surveys like LSST and Euclid. We present a method to infer the SEDs and spatial distribution of both the galaxy and its interstellar dust from LSST-like multi-band imaging. To stabilize this massively underconstrained inverse problem, we utilize two score-matching models as data-driven priors: the first informs our inference of the galaxy's underlying shape, the second informs the galaxy's dust morphology conditioned on the current estimate of the galaxy's shape. We believe that this is the first time a set of coupled score-matching models have been utilized to solve a complex data-analysis challenge. We demonstrate with realistic simulations that we can accurately measure the parameters of the underlying host and its dust content. In addition to providing galaxy SEDs unbiased by dust attenuation for subsequent analyses, such as photometric redshifts, our dust maps will allow the study of the interplay between star formation and dust production and destruction.